/  AD-A149  827 BOUNDHRV  ELEMENTS  FOR  OEBOND  STRESS  ANALYSIS;  FRACTURE  1/1 
MECHANICS  STRESS  A.  .  <U>  IMPERIAL  COLL  OF  SCIENCE  AND 
TECHNOLOGV  LONDON  (ENGLAND)  DEPT.  .  C  ATKINSON 
UNCLASSIFIED  26  APR  83  AFOSR-TR-84-1207  AFOSR-82-0121  F/6  12/1  NL 


AD- A 14 


RFPORT  DOCUMENTATION  PAGE 


’  ^  L  "  s.MBLt. 


AFOSR-TR-  3  4  .  1  207 


■  Z  GOVT  *:  CESSION  NO  3  CFCiPiLn^’S  CATALOG  M«>‘  R 


BOUNDARY  ELEMENTS  FOR  DEBOND  STRESS 
ANALY  SIS;  Fracture  mechanics  stress  analysis 


5  ’yp£  OP  PCPCRT  &  PERiDD  CTuP-fr 

FINAL 

01  MAR  82-28  FEB  83 

6  PERFORMING  0  RG.  REP' PT  Nc-MBE» 


7  AuThOR;*.- 


COLIN  ATKINSON 


contract  CP  GRAN*  n.MBlR  s 


AFOSR-S2-G131 


performing  organization  name  and  address 


ID-  PROGRAM  ELEMENT  PROJEC’ 
AREA  &  WORK  UNIT  NUM9ERS 


IMPERIAL  COLLEGE  OF  SCIENCE  k  TECHNOLOGY  61102F 
LONDON  SW7  2BZ  :  UNITED  KINGDOM  2  30  7/B2 


CONTROLLING  OFFICE  NAME  AND  ADDRESS 


12.  REPORT  DATE 

EOARD  26th  APRIL  19S3 

Box  14  13  NUMBER  or  PAGES 

FPO  NY  09510  _ _ _ _ _ SI 

I  MOMTGPiN  G  AGENCY  NAME  &  ADDRESS':/  different  from  Controlling  Office)  1  15  SECURITY  Ck_ASS.  (of  :h 

vir  Force  Office  of  Scientific  Research/NA  UNCLA8RIFI 

Soiling  AFB  DC  20332-6600  _  UU  1 

IS*  DECu  A5S1P1C  ATlON  D! 

schedule 


Dl  ST  Rl  9  u  T]  ON  STATEMENT  (of  ibis  Report, 


Approved  for  public  release;  distribution  unlimited 


KEY  WORDS  (Continue  on  reverse  side  it  necessary  ond  identity  by  block  number) 

Boundary  integral  equations 

BOUNDARY  ELEMENTS 
BOUNDARY  ELEMENT  METHODS 
FRACTURE  MECHANICS 
STRESS  SINGULARITIES.. 


20.  ABSTRACT  (Continue  on  reverse  side  It  necessary  and  Identify  by  block  number) 

To  understand  and  quantify  debonding  events  in  composite  rr.aterials 
’.t  is  necessary  to  have  accurate  stress  analysis  in  the  area  of 
stress  concentrators.  Boundary  integral  equation  (BIE)  techniques 
ire  proposed  as  a  potentially  more  accurate  approach  than  finite 
,  ilements.  In  particular,  the  debond  stress  analysis  of  a  rod  being 
pulled  out  of  a  matrix  material  is  considered.  A  model  problem  is 
posed  consisting  of  an  edge  cracked  bar  under  longitudinal  shear 
deformation  and  a  BIE  formulation  is  propused  using  a  subtraction 


DD  1  j  AN  ^TJ  1473  EDITION  OF  I  NOV  ES  IS  OBSOLETE 


UNCLASSIFIED 

SfCjHi’v  Classification  *  -  ?•  r  *-E 


058 


(Conjj  _ 


of  singularity  technique.  An  extensive  review  of  integral  eqvation 
approaches  to  various  fracture  mechanics  problems  is  also  prerented 
with  suggested  new  solution  techniques  for  unsymmetric  crack  confia 
urations.  ,, v 


/  1  f  I  u 


Accep'iion  For _ _ 

ktt:.  ^ 


.  co*  * 

mAPKi:  r* 

v  y 


♦  -  Codes 


SFOSR-TR-  i 4  _ 1 20  7 


Contract/Grant  No:  AFOSR  -  S2  -  0131 


BOUNDARY  ELEMENTS  FOR  DEBOND  STRESS  ANALYSIS: 
Fracture  Mechanics  Stress  Analysis 


Colin  Atkinson 

Department  of  Mathematics 

Imperial  College  of  Science  &  Technology 

LONDON  SW7  2BZ 


26th  April  1983 


Final  Report,  1st  March  1982  -  28th  February  1983 


Approved  for  public  release; 
distribution  unlimited 


Prepared  for:  UNITED  STATES  AIR  FORCE 

Air  Force  Office  of  Scientific  Research 
Building  410,  Bolling  AFB,  D.C.  20332 


and  EUROPEAN  OFFICE  OF  AEROSPACE  RESEARCH 

AND  DEVELOPMENT 
London,  England. 


85  01  14  058 


aih  mens  ovnr*  <r 

NOTICE  f.'F  r\-.' 

Thla  i.  ' 


_  -|  apr-  1 
Di.jf 


S!  w  •  •  w*l*l« 


INTRODUCTION 


In  order  to  understand  and  quantify  debonding  events 
in  composite  materials  under  stress  it  is  necessary  to  have  accurate 
stress  analysis  in  the  vicinity  of  stress  concentrators.  In  recent 
work  [1]  the  debond  stress  analysis  of  a  rod  being  pulled  out  of 
a  Solithane  matrix  was  considered.  The  methods  used  for  the  stress 
analysis  were  (i)  finite  element  calculations  of  the  energy  release 
rate  associated  with  the  debonding  process  and  (ii)  special 
analytical  solutions  valid  for  short  crack  initiation  at  key 
geometrical  features  of  the  rod-matrix  geometry.  In  this  second 
approach  (which  is  potentially  the  most  accurate)  a  key  inqredient 
is  the  stress  field  at  geometrical  discontinuities  lying  on  the 
rod-matrix  interface.  The  finite  element  method  yields  very  inaccurate 
results  at  such  discontinuities  unless  special  elements  are  used. 

A  popular  alternative  to  the  finite  element  method  is  the 
Boundary  Integral  Equation  (B.I.E)  method  but  it  too  may  require 
special  attention  to  the  stress  singularities  at  geometrical  dis¬ 
continuities.  In  the  papers  [2]  and  [3]  a  comparison  of  various 
methods  of  dealing  with  crack  tip  stress  singularities  was  made  for 
the  solution  of  a  model  problem  by  the  boundary  integral  equation 
method.  The  main  purpose  of  the  project  reported  here  was  to  extend 
this  work  to  other  stress  singularities  by  considering  the  model 
problem  discussed  below.  A  secondary  objective  was  to  attempt  to 
develop  a  theory  to  extend  the  BIE  applicability  to  unsymmetric  crack 
configurations.  Some  progress  on  this  latter  objective  has  been 
made  and  a  possible  method  is  suggested  in  the  Appendix  where  an 
extensive  review  of  integral  equation  approaches  to  various  fracture 
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mechanics  problems  is  also  presented.  In  our  original  proposal 
we  suggested  the  model  problem  discussed  below  noting  that 
experience  gained  from  treating  this  simpler  problem  is  potentially 
readily  transferable  to  the  elastic  situation. 


Model  Problem  (longitudinal  shear  deformation) 


The  normal  stress  o  is  assumed  zero  on  all  boundaries  except  the 

n2  K 

top  and  bottom  of  the  slab  (Fig.1)  for  which  we  assume  fixed  displace¬ 
ment  loading 


“3  =  u30  »  =  1 

u3  =  -°30  »  =  -1 


(2.1) 


Because  of  the  longitudinal  shear  (anti-plane  strain)  assumption  the 
only  non-zero  stresses  are 


(2.2) 


and  the  associated  equilibrium  equation  reduces  to 


Before  considering  the  BIE  implementation  of  the  problem  it  is 
useful  to  investigate  the  kind  of  singularity  which  exists  at  notch 
tip  A  in  terms  of  the  internal  angle,  2a,  of  the  notch.  Expanding 
the  displacement  field  in  terms  of  a  local  co-ordinate  system  (r,$) 
based  on  the  notch  tip  A  gives 

Oj  **  t\a.j  cosAS  -*■  B^  sinX$)  (; 

as  the  dominant  term  of  the  notch  tip  0  <  X  <  1.  Note  that  this 
expression  for  u^  is  an  exact  solution  of  the  field  equation  above, 
although  it  does  not  have  the  capability,  of  course,  of  satisfying 
all  the  boundary  conditions  of  the  problem.  Nevertheless,  the 
conditions  on  the  notch  sides  are  satisfied  exactly  if 

1  ^u3 

— —  —  r  0  on  ■&  -  ct ,  2tt  —  a  ( 

and  hence  we  require 

-A^  sinXa  +  B^  cosXa  =  0 
-A^  sin  A(2m-a)  +  B^  cosX(2iT-a)  =  0 


which  is  possible  provided 


sinXa  cosX(2n-a)  =  sinX(2ir-a)cosXa 


(2.7 


i.e.  sinX(2m-2a)  =  0. 

Thus  X(n-a)=  n/2  gives  the  dominant  singularity.  It  is  precisely 
singularities  of  this  more  general  kind  which  occur  at  geometric 
discontinuities  of  rigid  fibres  in  composite  materials. 

from  the  above  analysis  we  expect  that  in  the  vicinity  of  the 
notch  tip  (A  of  Fig.1)  the  dominant  part  of  the  displacement  u^  takes 
the  form 


n  X  /  TT 

u3  =  BQr  cos(^ 


($-ak 


(2.e; 


where  X  =  7t/2(tt-ol), 


(2.9 


and  the  constant  is  to  be  determined  from  a  full  numerical  analysis 
of  the  problem.  It  is  worth  noting,  for  later  use,  that 
expression  (2.8)  with  X  replaced  by  -X  also  satisfies  equation  (2.3) 
and  the  boundary  conditions  (2.5).  Although  the  behaviour  of  this 
solution  near  r  =  0  gives  infinite  values  for  u^  so  the  corresponding 
solution  is  non-physical,  it  is  nevertheless  useful  in  deriving  the 
physical  solution  as  we  shall  see  later. 


3. 


The  boundary  integral  equation  method 

The  starting  point  is  the  boundary  integral  equation 


f  T(X,Y)[u(Y)-u(X)lds  =  f  U(X,Y) 


3u(Y). 

—  as 


VX 


(3.i: 


*\  k 
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U(X,Y)  =  (  2?)  ln(1/r)  with  r  =  |V-X|  (3.2) 

2 

98  is  the  boundary  of  fl  a  bounded  open  region  in  R  . 

T(X,Y)  =  Vy  U(X,Y).N(V)  (3.3) 

where  N(_Y)  is  the  outward  unit  normal  at  the  point  ]f£38. 

In  equation  (3.1) 

J  f(X,Y)ds  (3.4) 

an  n  38(X.,€) 

and  the  region  fi(X^€)  is  similar  to  8  but  excluding  a  region  of  radius 
€  centred  on  the  point  X.  The  function  u  satisfies 


V2u 

It 

o 

in 

8 

u  = 

u 

on 

3fi1 

(3.5) 

3u 

3N 

l 

on 

3^2 

with  38  =  38  ^  u  • 

The  problem  is  thus  to  compute  approximations  u  and  t  to  u  and  t  on 
3^2  ar>d  38  respectively  from  the  integral  equation  (3.1)  and  boundary 
conditions  (3.3). 

Any  particular  boundary  element  method  is  characterised  by: 
(a)  the  choice  of  the  classes  of  functions  used  to  approximate 

u  and  t  and  (b)  the  way  in  which  the  integral  equation  (3.1)  is  used 
in  the  determination  of  u  and  t.  Approximating  u  and  t  (  =  3u/3N)  by 
linear  combinations  of  computationally  convenient  functions  v^u\ 


■f  f(X,Y)ds  =  lim 

J  -  y  €*o 
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(t) 


Z  =  1,...,N^and  v^,  ,  £'  =  1 , . . . , respectively  gives 


u  =  I  vi<u)  -  1  =  [ 


tZ'VZ 


(t) 


(3.6) 


JU 1 


Jl '  =1 


This  is  discussed  more  fully  in  section  4  of  the  appendix  where  an 
application  to  a  cracked  body  is  considered. 

To  extend  these  methods  to  notches  such  as  that  shown  in 
Fig.1  we  have  considered  the  following  three  methods. 


Method  1 .  Attempt  to  position  nodes  on  the  element  closest  to  notch 
tip  A  so  as  to  simulate  the  correct  displacement  behaviour  on  the 
notch.  From  equation  (2.4)  it  follows  that  the  required  displacement 
behaviour  is 

A  ,  ,  TT 

u  ~  r  where  A  =  - <r 

2 ( 7T  -  a) 


when  a  =  0,  the  notch  becomes  a  crack,  A  =  £  and  this  case  is 
considered  in  section  4.2.3  of  the  appendix  by  the  use  of  a  quarter 
point  node  and  quadratic  elements.  Our  initial  investigations  of 
the  case  \  i  \  suggest  that  this  approach  will  work  only  for  special 
values  of  A  e.g.  if  A  =  1/3  cubic  elements  would  probably  be  necessary. 

We  have  not  therefore  pursued  this  avenue  further  but  note  that  some 
work  on  this  approach  to  finite  elements  exists  see  (e.g.  Wait  [4]). 

Method  2.  Subtraction  of  singularity  technique.  Extending  the  argument 
given  in  section  2  in  deriving  the  dominant  term  for  u^  near  the  notich 
tip  one  can  write 


(3.7) 


u,  =  )  B .  r  cos(A .  ($-a) ) 

3  L  i  i 

i=0 

,  .  (2i+1)7i 

where  A .  =  —=7 — -~r 

1  2(7T-a) 

In  this  method  expression  (3.7)  is  subtracted  from  the  problem  before 
the  B.I.E.  solution  procedure  is  implemented.  Additional  equations 
for  the  unknown  coefficients  B^  can  be  determined  by  imposing  the  boundary 
condition  on  u^  at  additional  points  on  AC  (see  Fig.1).  On  AC  the 
boundary  conditions  u-j  =  0  follows  from  symmetry,  the  numerical 
solution  of  the  integral  equation  being  effected  on  the  boundary  ACDEF. 

This  method  has  been  used  effectively  for  a  crack  problem  discussed 
in  the  Appendix.  The  problem  of  a  re-entrant  corner  with  interior 
angle  3tt/2  has  been  calculated  by  this  method  and  good  agreement 
obtained  with  results  obtained  by  numerical  conformal  mapping  (Papamichael 
and  Whiteman  [5]).  A  method  similar  to  that  suggested  here  has  also 
been  considered  by  Symm  [6].  At  present  we  feel  that  this  method  is 
potentially  the  most  accurate  for  dealing  with  stress  singularities 
but  it  is  useful  to  calculate  the  coefficient  of  the  leading  term 
of  (3.7)  by  an  alternative  method  as  a  means  of  comparison. 

Method  3.  Evaluation  of  sharp  notch  stress  concentration  factors  by 
means  of  a  reciprocal  theorem. 

In  section  2  of  the  Appendix  methods  are  discussed  for  the 
evaluation  of  stress  intensity  factors  by  the  use  of  certain  invariant 
integrals  evaluated  far  from  the  crack  tip.  The  only  one  of  these  methods 
that  is  applicable  to  the  sharp  notch  problem  discussed  in  section  2 
when  a  i  0  is  that  based  on  Betti '  s  reciprocal  theorem  given  in 
section  2.2  of  the  appendix. 
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For  the  model  problem  of  section  2  this  can  be  written 

J  wu  ds  =  J  is  ds  (5-8) 

s  s 

where  u  =  and  u  satisfy  equation  (2.3)  and  S  is  envisaged  as  a 
closed  contour  such  as  shown  dotted  in  Fig.1.  The  divergence  theorem 
and  equation  ''2.3)  can  be  used  to  prove  equation  (3.8). 

If  for  tj  we  choose  the  solution 

‘  "X0 

u  =  r  cos  (X_  (■&-<! ) )  (3.9) 

with  A  =  ^  ,  ,  then  the  integrals  in  equation  (3.8)  are  zero 

when  taken  along  the  notch  sides  Af  and  AF^  because  of  the  boundary 

conditions  8u/3N  =  0  on  them.  The  integral  around  an  infinitesimal 

loop  enclosing  the  notch  tip  is  proportional  to  the  unknown  coefficient 

Bn  on  account  of  the  dominant  behaviour  of  the  displacement  field  u  =  u, 
U  X  ^ 

Ao 

there  (i.e.  u,~Bnr  cos(An($-a) ) ,  for  small  r  see  equation  (2.8)). 
j  u  U 

Thus  the  coefficient  can  be  determined  from  the  integrals  (3.8) 
evaluated  on  a  contour  far  from  the  notch  tip  e.g.  on  the  outer 
boundary  FEDC  D^E^F^,  the  values  of  u  and  9u/3NI  being  evaluated  there  by 
a  B.I.E.  method. 

A .  Concluding  remarks 

We  believe  that  method  2  discussed  in  section  3  is  the  most 
effective  way  of  computing  the  stresses  near  a  sharp  notch  tip  if  simple 
elements  are  used  to  model  the  boundary.  If  higher  order  elements  are 
used  more  accurate  results  could  presumably  be  obtained  for  special 
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cases  with  accurate  element  modelling  of  the  displacement  field  at  the 
notch.  More  work  is  required  to  make  an  extensive  comparison  of  all 
possibilities,  however,  at  present  we  recommend  method  2  together 
with  use  of  the  reciprocal  theorem  as  outlined  in  method  3. 


s 


References 


Atkinson,  C.,  Avila,  J.,  Betz,  E.  and  Smelser,  R.E.  (1982). 

The  rod  pull  out  problem,  theory  and  experiment,  J. 

Mech.  Phys.  Solids  30.»  97-120. 

Xanthis,  L.S.,  Bernal,  M.J.M.  and  Atkinson,  C.  (1981).  The 

treatment  of  singularities  in  the  calculation  of  stress 
intensity  factors  using  the  boundary  integral  equation 
method,  Comp.  Methods  in  Appl.  Mechs.  and  Engng., 

26,  285-304. 

Atkinson,  C.,  Xanthis,  L.S.  and  Bernal,  M.J.M.  (1981). 

Boundary  integral  equation  crack  tip  analysis  and 
applications  to  elastic  media  with  spatially  varying 
elastic  properties,  Comp.  Methods  in  Appl.  Mechs.  and 
Engng.,  29_,  35-49. 

Wait,  R.  (1977).  Singular  isoparametric  finite  elements, 

J.  Inst.  Maths.  Applies.  _20,  133-141. 

Papamichael,  N.  and  Whiteman,  J.R.  (1972).  A  numerical 

conformal  transformation  method  for  harmonic  mixed  boundary 
value  problems.  Brunei  University,  Dept,  of 
Mathematics,  Report  No.  TR/18. 

Symm,  G.T.  (1973).  Treatment  of  singularities  in  the  solution 


of  Laplace's  equation  by  an  integral  equation  method. 
NPL  Report.  NAC  31. 


FIGURE  . 


T 


APPENDIX 

1.  Introduction 

2.  Stress  intensity  factors  and  invariant  integrals 

3.  Integral  equation  methods  for  crack  tip  stress  analysis 

A.  Modelling  and  numerical  results 


5. 


References 


Vi 


1 .  Introduction 

In  this  report  the  boundary  integral  equation  (B.I.E.)  method  is 
discussed  with  particular  reference  to  its  application  to  fracture 
mechanics  stress  analysis.  Although  the  B.I.E.  method  has  been  us<d  to 
treat  a  variety  of  boundary  value  problems,  those  encountered  in  fracture 
mechanics  possess  certain  unique  features  which  require  special  treatment. 
There  are  essentially  two  distinct  difficulties  which  are  encountered 
when  applying  standard  B.I.E.  procedures  to  fracture  models  in  which  a 
crack  is  modelled  as  having  a  planform  of  zero  thickness  (e.g.  a  line  crack 
in  two  dimensions). 

These  are:- 

(i)  The  indeterminacy  encountered  when  (in  plane  strain  or  stress),  the 
direct  B.I.E.  formulation  of  a  thin  ellipse  is  used  to  represent  a 
crack  when  the  ellipse  degenerates  to  a  line.  A  similar  feature 
occurs  in  three  dimensional  problems.  We  describe  this  degeneracy 

more  fully  later  in  this  introduction,  attempts  to  remove  this  difficulty 
are  described  in  section  3. 

(ii)  The  existence  of  a  stress  singularity  at  a  sharp  crack  tip  which, 
requires  accurate  boundary  element  modelling  in  order  to  obtain 
reliable  numerical  results.  Such  modelling  is  described  in  detail 
in  section  4.  Note  however  that,  although  a  square  root  stress 
singularity  occurs  in  the  elastic  analysis  of  a  crack  in  a  homogeneous 
medium,  other  stress  singularities  may  be  encountered  (see  [1)  for  a 
review) . 

It  is  worth  noting  that  information  about  ’ s t res s- i n tens i ty  far  tots' 
(notation  to  be  explained  more  fully  in  section  2)  can  be  obtained  frrr. 
stress  and  displacement  fields  far  from  the  crack  tip  by  judicious  ucr>  of 
certain  path  independent  integrals.  Such  integrals  are  discussed  in 
section  2  and  their  use  will  be  illustrated  in  section  4. 
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To  understand  the  first  difficulty  discussed  above  we  review  here 
the  derivation  of  the  direct  integral  formulation  of  Rizzo  [2] .  The 
boundary  integral  equations  can  be  deduced  via  an  appropriate  reciprocal 
theorem,  fundamental  solution  and  certain  limiting  processes  as  field 
points  approach  the  boundary.  For  an  elastic  medium  the  starting  point  is 
Betti's  second  theorem  which  states  that  for  continuous,  finite  stresses 
and  zero  body  force 


f  t.u.*dS  =  f  t.*u.c 

1  i  l  J  l  l 


(8.1) 


where  t^  *  o„  n^  is  the  traction  vector,  u^  is  the  displacement  vector 
and  these  together  with  the  field  with  superscript  *  denote  equilibrium 
stress  states  which  are  well  defined  at  all  points  in  a  region  R  with 
bounding  surface  S.  In  the  above  reciprocal  theorem  the  stress  state  with 
superscript  *  will  be  chosen  to  correspond  with  the  fundamental  solution  to 
Navier's  equations  which  for  isotropic  elasticity  can  be  written:- 


Fundamental  solutions  for  two  dimensional  plane  strain  ( u  is  the  shear 
modulus,  v  Poisson's  ratio) 


u •  * ( P ,q)  *  [(3-4v)log(l/r)6.  .  +  r  .r  ,]e  ./8ttu(1-v) 
1  1 J  » 1  >  J  J 


(8.2) 


t.*(p,q) 


-  —  [(l-2v)6..  +  2r  .r  .]  -  (l-2v)[n.r  .-n.r  .  ]  }e  . /4ir(  1- v) 

r  dn  ij  ,i  ,j  J  ,i  i  ,J  j 


(8.3) 


where  r(p,q)  is  the  distance  between  points  p(x)  with  co-ordinates  x.  and 

of  q(x)  with  co-ordinates  y-,  then  r  .  is  defined  by 

1,1 


(8.4) 


the  e.  are  unit  vectors.  The  tractions  on  an  arbitrarv  surface  su  r  t  mind  i  nj-, 

J 

the  point  p(x)  can  be  computed  from  u.*(p,q)  by  differentiation  of  p(x) 
and  are  given  in  equation  (8.3)  above  where  the  normal  is  taken  at  q(x). 

Apart  from  being  a  solution  of  the  equilibrium  equation  when  p  f  q  the 
above  solution  has  the  significant  property  that  on  a  small  disc  Sr  of 
radius  e  centred  at  the  point  p(x) 

lim  /  t.*dS  =  6..e.  .  (8.3) 

e-0  1  1J  J 

Similarly  in  three  dimensions  we  have: 

Fundamental  solution  in  3  dimensions 

u.*(p,q)  =  — - r-  [(3-4v)6..  +  r  .r  .  ]e  ./16nu(l-v)  (8.6) 

l  r(p,q)  ij  ,i  ,j  j 

and 

t .  *  ( p ,  q )  =  -  -4r  {^[(l-2v)6..  +  3r  .r  .]  -  (l-2v)(n.r  -n.r  .  )  )e  . /8-i  ( 1  -  ) 
1  r2  dn  lj  ,i  .J  j  i  ,j  l 

(8.71 

where  e.  are  unit  vectors  in  the  co-ordinate  directions. 

J 

A  significant  property  of  this  solution  is  that  on  a  small  sphere 
of  radius  e  about  the  point  p(x) 


lim 

e-0 


J 


S 

t 


t . *  dS  *  5 .  .e  . 
i  U  J 


(8.8) 


Thus  the  fundamental  solution  corresponds  to  the  well-known  Kelvin 
problem  of  a  point  load  in  an  infinite  body.  Using  the  above  one  can  write 


u.*(p,q)  *  U..(pfq)e.;  t.*(p,q)=T..(p,q)e.  (8.°) 

l  ji  j  l  Ji  J 

such  that  the  first  index  in  U,.(p,q)  and  T^(p,q)  corresponds  to  the 
direction  of  the  point  load  and  the  second  index  refers  to  the  component 


of  the  respective  displacements  and  tractions. 


4 


Since  the  above  fundamental  solutions  are  singular  at  p(x)  we 
apply  equation  (8.1)  to  the  region  R  with  a  region  R  (a  ball  of  radius  e 


around  p(x))  deleted  from  it.  If  this  ball  has  surface  equation  (8.1) 


is  rewritten 


/  t .u. *dS  =  /  u . t .  * 

'  „  l  l  ii 


dS 


(8.10) 


S+S 


S+S 


The  expressions  (8.3)  and  (8.4)  or  (8.6)  and  (8.7)  for  u  .*  and  t* 


are  now  substituted  into  equation  (8.10)  and  the  stress  field  o..(p) 


and  displacement  u^(p)  assumed  continuous  and  bounded  for  any  point  p(x) 
belonging  to  R.  The  following  limiting  results  are  also  deduced 


lim 

e+0 


t.u.*dS  = 
i  l 


(8.11) 


and 


lim  /  u.t.*dS  =  u.(p)s..e.  =  u.(p)e; 

e-*>  s  11  1  *  >  1 

e 


(8.12) 


Using  these  expressions  in  equation  (8.10)  and  taking  the  limit  z  •+  0 
leads  to  a  standard  result  at  an  interior  point  p 


u  .  ( p )  =  J  t.(Q)U..(p,Q)dS(Q)  -  /  u.(Q)T..(p,Q)dS(Q) 

J  s  1  J1  s  1  J1 


(8.13) 


where  the  tensors  U..  and  T..  can  be  deduced  from  the  definitions  (8.9). 

Ji  Ji 


The  above  equation  thus  gives  a  representation  for  the  displacements  at 
any  interior  point  p(x)  and  the  stresses  may  be  calculated  by  differentiation 


with  respect  to  the  field  point  p(x)  with  co-ordinates  x^ .  The  result  can 


be  written 


o.jtp)  -  /  tk(Q)Dk..(p,q)dS(Q>  -  /  uk(Q)Sk..(p,Q)dS(Q) 

s  s 


(P.14) 


*"*  •"*  j  «*•  ,* 


1 

1 

i 

i 


where  the  tensors  D,  .  .  and  oi,  .  .  art  determined  from  U..  and  T.. 

kij  kxj  )i  ji 

respectively,  by  differentiation  with  respect  to  the  field  point  p(x)  and 
using  the  stress-strain  relations. 

The  above  results  give  the  stress  and  displacement  at  an  interior 
point.  The  well  known  boundary  integral  equation  is  obtained  by  taking 
the  limiting  form  of  equation  (8.13)  (sometimes  known  as  Somigliana's 
identity  for  the  displacement  vector)  and  letting  the  interior  point  p(x) 
tend  to  a  boundary  point  P(x).  Taking  due  account  of  the  discontinuity 
of  the  second  integral  on  the  right  hand  side  of  equation  (8.13)  as 
p($)  -*■  P(x)  from  the  interior  leads  to  the  integral  equation 

Uj(P)/2  +  /  u.(Q)T. .(P,Q)dS(Q)  =  /  t.(Q)U. . (P,Q)dS(Q)  (8.1") 

S  S 

The  difficulty  with  the  application  of  the  standard  B.I.E.  formulation 
to  fracture  mechanics  has  been  outlined  clearly  by  Cruse  [3,4,5],  we 
follow  his  derivation  here.  He  considers  the  two  dimensional  situation 
shown  in  figure  1  where  S,  P  and  f  represent  the  regular,  upper  and  lower 
crack  surfaces  as  shown.  The  Somiglinna  identity  (equation  (8.13))  lor 
the  displacement  at  an  interior  point  p(x)  is  thus  given  by 

u  .(p)  =  J  U . . (p,Q)t . (Q)dS(Q)  J  T  .  . (p,Q)u. (Q)dS(Q)  (8.16) 

J  s+r  +r"  JX  1  s+r++r  J1  1 

The  kernels  U..  and  T..  are  defined  in  equation  (8.9)  and  have  the 
jl  jx 

property  that 

T-Cp.Q*)  -  -T.^(p,q")  (8.17) 

and 

U.  . ( p ,Q+)  «  +U.  .(p,q”)  (8  |S) 

lj  r  lj 

for  Q  and  q  lying  on  the  crack  faces  P+  and  T  respectively  except  neat 
the  small  crack  tip  region.  Note  that  the  sign  change  in  (8.1"',  i «  rim- 
to  the  opposite  normal  directions  of  the  two  crack  surfaces. 


•V*  A  «*"  »4’  .*•  .‘v*  .*• 


1  «’  t  t  i  ue  ll 


two  crack  surfaces  collapse  to  the  plane  T  (with  the  same  normal  as  1+) 
leads  to  the  equation 


u.(p)  =  /  U..(p,Q)t  (Q)dS(Q)  -  J  T..(p,Q)u.(Q)dS(Q) 

J  *  C  J  ^  * 


-  /  T^(p,Q)Au^(Q)JS(Q)  ♦  /  IK.(p.Q)  ItAQ)dS(Q) 


(8.19) 


where 


Au.(Q)  =  u.(Q+)  -  u.(Q~)  (8.20) 

J  J  J 

Zt j (Q)  =  t j(Q+)  +  t^(Q')  (8.21) 

»  I  *■»  •  •  *  j 

In  general  situations  are  encountered  where  the  crack  is  loaded  by  equal 
and  opposite  tractions  such  that  It.(Q)  =  0.  The  usual  next  step  is  to  let 
p(x)  •+•  P(x)  a  boundary  point  in  equation  (8.19)  in  order  to  obtain  the 
boundary  integral  equation.  However,  when  P(x)  is  a  point  on  the  plane  T 
(not  at  the  edge  of  T )  the  equation  becomes 


u.(P)  -  Au . (P)/2  +  /  T. . (P,Q)u. (Q)dS(Q) 

J  J  s  i 

♦  j  T. ,(P,Q)Au.(Q)dS(Q)  =  /  U..(P,Q)t.(Q)dS(Q) 
p  J1  i  s  J1  1 


(8.22) 


where  the  usual  term  for  the  jump  in  the  boundary  integral  has  been  taken 

and  the  integral  over  T  is  a  principal  value  integral.  Equation  (8.22) 

is  deficient  in  two  critical  respects.  In  the  first  place  while  a  single 

surface  P  is  being  treated,  two  variables  are  unknown.  These  are  Au^ 

and  Zu .  *  u.(P)  -  Au.(P)/2.  The  second  deficiency  is  that  equation  (8.22) 
J  J  J 

is  the  same  for  any  set  of  equal  and  opposite  crack  boundary  tractions 
(i.e.  such  that  Zt^(Q)  *  0  where  Q  lies  on  F ) .  Thus  this  integral 
equation  cannot  deal  with  situations  where  the  crack  is  loaded  internal ly. 

As  a  consequence  of  the  above  mentioned  deficiencies  inherent  in  a 
straightforward  application  of  the  boundary  integral  equation  (8.22)  to 


flat  crack  modelling,  a  variety  of  methods  have  been  suggested  which  each 


deal  with  special  cases.  No  all  purpose  method  applicable  to  a  variety 
of  unsymmetric  crack  plan  forms  seems  to  be  available  at  the  present  time. 
However,  a  possibility  is  suggested  in  section  3,  but  a  practical  test  of 
its  applicability  is  yet  to  be  made.  A  discussion  of  the  variety  of 
methods  which  have  been  suggested  for  special  situations  will  also  be  made 
in  Section  3.  Note,  however,  that  the  above-mentioned  difficulties  are 
all  specific  to  the  sharp  crack  model  of  fracture  mechanics.  Useful 
numerical  results  might  still  be  obtained  by  modelling  the  crack  as  an 
elliptical  hole  to  give  an  accurate  stress  field  far  from  the  crack  and 
then  use  invariant  integrals  to  calculate  the  (sharp)  near  crack  tip 
stress  field.  A  discussion  of  such  invariant  integrals  is  made  in 


sections  2.2  and  2.3. 


Stress  intensity  factors  and  invariant  mteerals 


As  discussed  in  the  introduction  there  are  certain  features  of  the 
stress  anslysis  of  theoretical  fracture  models  that  distinguish  these* 
problems  from  others  to  which  B.I.E.  methods  have  been  applied.  For 
completeness  we  include  here  a  brief  description  of  some  of  the  most 
commonly  used  notations  and  results,  notably  those  pertaining  to  a 
crack  in  a  homogeneous  elastic  body.  For  situations  where  cracks  meet 
interfaces  between  different  elastic  media  or  lie  along  a  bimaxerial 
interface  different  singularities  are  encountered,  for  a  review  of 
some  of  these  situations  see  [1)  . 


2.1  Fracture  modes  and  stress  intensity  factors 


A  key  feature  m  the  stress  analysis  of  the  crack  models  considered 
here  is  the  singular  stress  field  at  the  crack  tip.  For  a  crack  in  a 
homogeneous  medium  an  eigenvalue  analysis  (see  for  example  Williams  [6] ) 
gives  a  result  of  the  form 


o  ^  KF(6)r  +  non-singular  terms 


(8.23) 


where  (r,e)  are  polar-co-ordinates  based  at  the  crack  tip  (see  Fig.  2), 

F(9)  is  a  function  determined  by  the  eigenvalue  analysis  and  is  independent 
of  the  applied  load  or  body  shape.  Thus,  close  to  the  crack  tip  the 
interaction  with  the  applied  stress  field  is  fully  characterised  by  the 
coefficient  K.  This  coefficient,  which  depends  on  the  applied  load,  the 
shape  of  the  body,  and  the  crack  length,  is  called  the  stress  intensity  factor 
and  is  a  cornerstone  of  contemporary  fracture  mechanics.  In  general,  when 
the  applied  stress  induces  stresses  or.  the  crack  which  are  a  combination 
of  pure  tension  and  shear,  there  are  three  stress  intensity  factors 
associated  with  the  different  modes  of  fracture.  These  are  illustrated  in 


Figs  3.  A  standard  notation  is  to  use  K ^  for  the  opening  (mode  I).  K 
for  in-plane  sliding  (mode  II),  K  for  anti-plane  sliding  (mode  111) 
modes  of  relative  displacement  of  the  crack  surfaces.  The  following 
results  are  standard  for  isotropic  bodies  in  plane  strain,  plane  stress 
or  anti-plane  strain.  For  the  pure  opening  mode,  we  have 


cos(8/2) 


( 2ur ) ' 


1  -  sin( 0/2 )sin(36/2 ) 1 
sin(6/2)cos(39/2) / 
1  +  sin( 0/2)sin(3e/2 )  | 


(8.24) 


for  the  singular  near  crack  tip  stress  field,  and 


„  cos(0/2 )  (<-1+2  sin  (0/2)] 

h  (_L)%  •  2 

2u  2n  !  sin(e/2)[<+l  -  2  cos  (6/2)] 


(8.25) 


where  <  =  3-4v  (plane  strain)  <  =  (3-v)/(l+v)  (plane  stress). 

Note  =  v(°h+022^  ^or  P^ane  strain.  For  the  in  plane  sliding  mode 


(2nr ) 


!-sin(6/2)(2  +  cos  (0/2)cos(  39/2)] 
cos(6/2)[l  -  sin( e/2)sin(30/2 ) ] 
sin(9/2)cos(e/2)cos(30/2) 


'6.26) 


|  1  l  =  —  (— )* 
k  ^  2” 


sin(0/2)[>r  +  1  +  2  cosz(0/2)] 
-cos(6/2)[<  -  1  -  2  sin2(0/2)] 


(8.27) 


For  anti-plane  sliding  we  have 


(2wr)1/2  )  cos( 6/2 ) 


-s in( 0/2  ) 


(«.?«) 
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and 

g 

u,  -  2  (y-) 1/2  sinCe/2)  (8.29) 

0  y  /  TT 

In  each  of  the  abo  \e  v  is  Poisson's  ratio,  y  is  the  shear  modulus , 

(r,0)  are  polar  co-ordinates  based  at  the  crack  tip  as  shown  in  Fig.  2, 
and  the  crack  lies  on  the  plane  Xj  =  0  as  shown  in  Fig.  2. 

It  is  obvious  from  the  above  formulae  that  the  stress  intensity 
factor  K  is  an  important  parameter  characterising  the  singular  near  crack 
tip  stress  field.  It  connects  with  fracture  theory  through  the  postulate 
that  fracture  will  commence  (i.e.  a  crack  will  grow)  when  the  stress 
intensity  factor  reaches  a  critical  value  Kc ,  this  value  being  an 
experimentally  determined  material  property  assumed  independent  of  the 
shape  of  the  body.  Such  a  fracture  criterion  is  however  unambiguous  only 
when  a  single  mode  is  operative.  In  a  mixed  mode  situation  which  occurs 
for  example  when  a  crack  is  aligned  at  an  angle  to  an  applied  tensile 
stress  there  is  still  discussion  as  to  what  is  the  correct  fracture 
criterion  to  use.  We  concern  ourselves  in  this  article  with  accurate 
procedures  for  calculating  the  near  crack  tip  stress  field  for  either  single 
or  mixed  mode  situations.  The  question  as  to  what  fracture  criterion  is 
best  suited  to  predicting  mixed  mode  fracture  will  not  be  addressed, 
however,  see  [1]  for  a  review  and  discussion  of  this  point. 

In  addition  to  the  above  "critical  stress  intensity  factor"  fracture 
criterion  there  is  the  'energy  release  rate'  fracture  theory  of  Griffith. 
This  states  that  for  a  pre-existing  crack  to  grow,  the  decrease  of  the 
total  energy  (elastic  energy  plus  potential  energy  of  the  loading  mcrh m i cm) 
must  be  equal  to  or  exceed  the  surface  energy  of  the  two  newly  created 
fracture  surfaces.  The  energy  release  rate,  or  crack  extension  for<°,  C, 
can  thus  be  defined  as  the  amount  of  energy  released  from  the  s\«tf'm 


i 


(cracked  specimen  plus  loading  mechanism)  for  unit  advance  of  the  crack 
front.  For  fracture  in  cither  of  the  three  modes  shown  in  Kip  J 
it  can  be  shown  that 


K^U-v2) 


Kn2n-,2> 


(8.30) 


where  E  is  Young's  modulus.  The  Griffith  fracture  criterion  then  says 
that  fracture  occurs  when  G  =  2y  (y  being  the  surface  energy  of  each  of 
the  newly  created  surfaces).  Thus  when  a  single  mode  is  operative  the 
criteria  of  critical  stress  intensity  factor  and  critical  energy  release 
rate  are  equivalent.  However,  in  the  mixed  mode  situation  for  a  crack 


growing  in  its  own  plane  the  energy  release  rate  is 


G  =  (K. 


.  K  2) 

II  E 


(6.31) 


and  no  simple  association  between  the  two  criteria  is  possible 
(cf  .  reference  [1] ) . 

The  above  discussion  applies  strictly  to  elastic  media  and  brittle 
solids  but  nevertheless  has  had  quite  widespread  practical  success. 

The  subject  of  elastic-plastic  fracture  mechanics  is  a  modification 
of  the  theory  described  above  designed  to  t3ke  into  account  plastic 
deformation  at  a  crack  tip.  A  detailed  description  of  this  subject  is 
outside  the  scope  of  the  present  article,  however  the  following  features 
are  noted  (see  [1]  or  [7]  for  a  more  complete  account.)  For  problems  of 
small  scale  yielding  the  singular  solutions  outlined  in  equations  (8.24) 
to  (8.29)  are  useful  in  setting  an  'outer'  approximation  to  the  stress 
field.  In  models  where  plastic  flow  occurs  by  slip  bands  emanating  from 
a  crack  tip  such  as  envisaged  by  Bilby  and  Swinden  (8)  or  by  suprrd i s 1 ot a t i 
modelling  of  plasticity,  Atkinson  and  Kay  [9]  (see  also  |10]),  it  is 
important  to  have  some  information  of  higher  order  terms  in  (8.231  in  nrder 


to  determine  the  shear  stress  induced  on  the  slip  band.  As  we  shall 
see  later  this  information  is  readily  available  with  certain  kinds  of 
element  modelling  in  the  B.I.E.  integral  equation. 

Another  concept  much  used  in  elastic-plastic  fracture  mechanics  is 
the  J-integral  introduced  by  Rice  (11).  This  integral  is  identical  to 
a  certain  integral  of  the  energy-momentum  tensor  discussed  by  Eshelby  (12), 
in  the  context  of  generalised  forces  on  point  defects  and  inhomogeneities 
in  elastic  fields.  Under  certain  circumstances  integrals  based  on  the 
energy  momentum  tensor  are  invariant  and  by  suitably  deforming  the  contour 
of  the  integral,  information  about  near  crack  tip  stress  fields  can  be 
deduced  from  numerical  information  obtained  on  the  boundary  of  a  region 
removed  from  the  crack  tip.  A  brief  derivation  of  some  of  these  integrals 
will  be  given  in  the  next  section  as  a  precursor  to  their  numerical 
implementation  in  section  4. 

2.2  Invariant  integrals  based  on  the  energy  momentum  tensor 

We  have  already  indicated  briefly  in  section  2.1  how  certain 
invariant  integrals  can  be  used  as  tools  for  deducing  numerical  values 
stress  intensity  factors.  Here  a  simple  deduction  of  integrals  which  can 
be  derived  from  a  (perhaps  pseudo)  energy  momentum  tensor  is  made.  The 
prescription  used  is  similar  to  that  described  by  Eshelby  [13)  although 
equivalent  results  can  be  obtained  by  a  formal  application  of  Noether's 
theorem  (see  Knowles  and  Sternberg  [14]  and  for  somewhat  equivalent  earlier 
work  by  Gunther  [15]). 

We  suppose  a  Lagrangian  (L)  is  given  which  has  the  property  that 
the  conditions  of  stationarity  of  the  functional  / L  dv  (where  v  is  the 
volume  of  the  body)  lead  to  Euler  equations  which  are  themselves  the  fi°ld 
equations  of  the  continuum  being  studied.  In  order  to  retain  some  generality 


suppose 


L  ■  L  ( X .  ,  u .  ,  u  .  .) 

i  i  i  .J  l  l  ,  J 


1 '» 


(K.32) 


thus  L  depends  on  position  (cartesian  co-ordinates  X^)  and  two  independent 

vector  fields  u. ,  $ .  and  their  gradients  u.  .  =  3u . / 3X .  etc.  .  The  Euler 
11  i ,  J  i  j 

equations  are  thus 


9  ,  9L  >  9L 

9X .  l9u.  .  9u.  ~  °  and  3X .  ' 3$ .  /  36. 

J  i,J  1  J  i.J  1 


3  ^  3L  ^  _  3L 


=  0 


(8.33) 


and  the  summation  convention  with  respect  to  repeated  indices  (i  =  1,2,3) 
has  been  used.  If  we  now  define  the  tensor 


Pj* 


9L  ,  9L 

-  u  +  -  <+>  —  Lc 

3ik  .  i,£  3$ .  .  yi  ,4  j  l 


(8.34) 
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then  it  can  be  shown  by  direct  calculation,  using  the  above  Euler  equations, 
that 

9P.„ 


j  ^  /  3L  \ 

9X .  “  "  3X„  exp 

J  5. 


(8.3S) 


where  (3L/3X  )  means  that  all  variables  are  held  constant  except 
£  exp 

explicit  dependence  on  X^.  This  last  equation  is  important  because  it 
means  that  if  we  define  the  integrals 


F  =  J  P.„  dS. 
*  S  J 


(8.36) 


where  S  is  a  specified  surface,  then  if  the  surface  integral  encloses  a 
volume  in  which  there  are  no  singularities ^use  of  the  divergence  theorem 
and  equation  (8.35)  reduce  the  corresponding  volume  integral  to  zero 
provided  L  does  notdepend  explicitly  on  X^ .  This  means  that  a  surface 
integral  taken  around  a  crack  tip  can  be  related  to  an  integral  far  from 
the  tip  and  since  under  favourable  circumstances  the  nesr  tip  integral 
can  be  calculated  directly  in  terms  of  a  single  unknown  stress  intensity 
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factor  evaluation.  It  should  be  stressed  however  that  this  is  only 

possible  provided  the  appropriate  components  of  P  are  zero  along  the 

j  £ 

crack  surface.  When  this  last  condition  is  satisfied  the  near  tip 
integral  will  be  exactly  equal  to  the  integral  evaluated  far  from  the 
crack  tip  if  L  does  not  depend  explicitly  on  the  spatial  co-ordinate  X  . 

We  illustrate  the  above  procedure  with  two  fairly  general  examples 
(the  standard  example  of  classical  elasticity  will  be  included  as  a 
special  case).  The  first  example  we  consider  is  that  of  micropolar 
elasticity  (cf  Atkinson  [16]  for  a  more  complete  account). 


2.2.1  Application  to  Micropolar  Elasticity. 

The  theory  of  micropolar  elasticity  is  reviewed  by  Eringen  [17]  and 
depends  on  the  usual  elastic  strain  tensor  together  with  an  independent 
'microrotation  vector' .  For  a  discussion  of  the  limitations  of  the  classical 
couple-stress  theory  and  its  relation  to  the  micropolar  theory  the  reader 
is  referred  to  [17] . 

We  suppose  there  exists  a  strain  energy  density  function 


W  =  W(u .  ,  ,  )  (8.37) 

i.k  j  j,s 

where  is  the  micro-rotation  vector  and  a  comma  denotes  partial 
differentiation  with  respect  to  the  cartesian  co-ordinates  (X^.Xj.X^). 

We  assume  that  W  has  the  property  that 
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Replacing  L  by  -W  in  equations  (8.33)  above 
for  minimising  the  functional  j W  dV.  In  t 


couple  stresses  (m^  )  defined  above,  these 
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leads  to  the  Euler  equation 
erms  of  the  stresses  (t  .  )  and 

P  r 

equilibrium  equations  can  be 
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written 


The  energy  momentum  tensor  defined  in  equation  (8.34)  is  thus 
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(8.40) 


where  we  note  that  the  first  index  denotes  the  normal  direction,  i.e.  tjf 
is  the  stress  in  direction  2  acting  on  a  plane  of  normal  j.  The  rest  of 
the  theory  discussed  above  also  applies  to  this  case.  For  example  for  a 
plane  crack  the  crack  extension  force  or  energy  release  rate  (G)  can  he 
defined  as 


G  =  F  =  /  P .  dS .  (8.41 ) 

1  S  J*  J 

where  S  is  a  surface  enclosing  the  crack  tip  (in  the  2d  case  this  is  a 
cylinder  with  generators  parallel  to  the  X^  axis  so  that  the  integral  in 
(8.41)  is  effectively  a  line  integral  in  the  (X^,X9)  plane).  Note,  also 
that  for  a  stress  free  crack  the  tractions  t^  and  the  couple  stresses  m? . 
will  be  zero  so  the  integrand  P,^  is  zero  on  the  crack  faces  and  the 
integral  defining  G  taken  around  the  crack  tip  for  an  edge  crack  can  he 
deformed  into  the  far  field  provided  only  that  the  energy-density  W 
does  not  depend  explicitly  on  Xf  (compare  equation  (8.33)). 

Thus  for  micropolar  elasticity  the  invariant  integral  F^  can  be 
used  to  determine  the  energy  release  from  numerical  information  derived  at 
points  remote  from  the  crack  tip.  It  should  be  noted  however  that  it  is  not 
possible  in  general  to  extract  the  crack  tip  stress  intensity  factor  f’-nm 
this  information  since  the  energy  release  rate  contains  a  couple  stress 
intensity  factor  contribution  in  addition  to  the  classical  elastic 
contribution.  For  more  information  on  this  point  and  a  detailed  analysis 
of  the  internal  crack  problem  in  micropolar  elastic  media  see  Atkinson  and 
Leppington  [18]. 


The  above  analysis  applies  also,  of  course,  to  classical  elasticity 
when  is  identically  zero,  the  stress  tensor  t  then  being  symmetric 
in  the  linear  theory. 

2.2.2.  Coupled,  time  dependent,  linear  thermoelasticit;.  and  thermovisco¬ 
elasticity 

Recently  Atkinson  and  Smelser  [19]  have  applied  a  procedure  similar 
to  that  outlined  above  to  the  equations  of  coupled  time  dependent  thermo¬ 
viscoelasticity  under  conditions  applicable  to  stationary  cracks  disturbing 
time  dependent  temperature  and  stress  fields.  Initial  conditions  are 
considered  in  which 


0(t)  =  u.(t)  =  o..(t)  =  0  for  t  <  0  (8.42) 

i  lj 


u.  and  o. .  are  the  usual  displacement  vector  and  stress  tensor,  6( t ) 
i  U 

denotes  the  infinitesimal  temperature  deviation  from  the  base  temperature 


T  .  The  formulation  described  in  [19]  begins  by  Lapl ace-trans forming  the 
equations  of  motion  etc.  which  become 
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e.  .  ■  i  (u.  .  ♦  u.  . )  (8.43) 
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o.  .  =  p  G.  „  e,  „  -  p  <4.  .  "5  (anisotropic) 
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and  for  the  temperature 
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where  the  Laplace  transform  is  defined  by 


f(p) 
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e  ^  f(t)dt  . 


In  general  for  viscoelastic  media  the  coefficients 
functions  of  p  the  transform  variable. 


(8.44) 


(8.45) 


etc .  will  he 


The  above  field  equations  can  he  generated  from  a  Lagrangian 
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defined  as 
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where  t .  .  =  p  C .  „e, 
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(8.47) 


A  'pseudo'  energy  momentum  tensor  is  then  defined  as 


P  .  =  U  +  e  -  L6  . 
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so  that  we  deduce  that 
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and  the  integrals 
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(11.50) 


follow  as  described  earlier.  A  useful  property  of  the  integral  in 
this  case  is  that  provided  either  the  temperature  0  is  constant  or  the 
flux  k.„0  .  is  zero  on  the  crack  faces  then  P,„  is  zero  on  the  crack  face 

i2  ,i  12 

for  a  stress  free  crack.  Thus  in  this  case  the  integral  can  be  deformed 
into  the  far  field  as  discussed  earlier.  Also  the  near  field  integral, 
a  small  contour  round  the  crack  tip  can  be  explicitly  evaluated  in  t'MiDi 
of  the  coefficients  of  singular  transformed  stresses  o„  and  temperature 
gradients  0  In  favourable  circumstances  explicit  determination  of  three 

coefficients  can  be  made,  see  [19]  for  some  applications  of  these  results. 


2.2.3  L  and  M  integrals  in  elastostatics 

In  the  above  we  have  discussed  the  integrals  associated  with  the 

tensor  P .  paying  particular  attention  to  the  integral  F  which  can  be 
ij  1 

shown  to  be  related  to  the  energy  release  rate  G.  Under  certain  conditions 

other  invariant  integrals  can  be  derived  from  P„.  and  are  sometimes  useful. 

Ij 

We  describe  these  briefly.  In  addition  to  the  F^  integrals  obtained  in 
section  2.2,  Cunther  [15]  obtained  for  el astostatics  the  path-independent 
integrals 


I  ‘Vij  *  Vlj>  dSj 


(8.51) 


and 
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(8.52) 


where  S.,j,k  take  the  values  1.  2  and  3  and  P  *  8W(u  )/3u  .  is  the 

J  m,n  i,j 

Boussinesq  or  first  Piola-Kirchof f  stress-tensor  which  gives  the  component 

parallel  to  the  (rectangular)  co-ordinate  axis  of  the  force  on  a  surface 

element  which  was  of  unit  area  and  perpendicular  to  the  axis  before 

deformation.  The  energy  momentum  tensor  is  defined  as  in  equation  (8.34) 

with  L  *  -W  and  the  energy  density  W  is  a  function  of  u  *  9u  /3X  with 

m,n  m  n 

n  *  1,  2  or  3,  cijjt  *8  t*ie  permutation  tensor. 

In  section  (2.2)  it  is  shown  that  the  integrals  F^  will  be  path 
independent  if  L  “  -W  does  not  depend  explicitly  on  X^..  The  deformation 
may  be  non-linear  and  the  material  may  be  non-linear  but  the  energy  density 


must  be  homogeneous  if  each  of  the  F^  is  to  be  path  independent.  For  the 
integrals  these  conditions  are  the  same  and  in  addition  the  material  must 


be  isotropic. 


In  a  two  dimensional  state  of  plane^(8.52)  reduces  to 


-  1  x.  P,.  dS. 

S  1  £j  J 


M 


(R.S3) 


taken  along  a  plane  curve  S  with  normal  (n  .n^).  Summation  over  >.  and  j 
is  now  only  over  1  and  2.  An  application  of  the  F  and  M  integrals  to 
the  case  of  anti-plane  strain  is  made  in  section  4  to  obtain  crack  tip 
stress  intensity  factors  from  numerical  information  far  from  the  crack  tip. 


2.3  Invariant  integrals  deduced  from  Betti's  reciprocal  theorem 


The  integrals  discussed  in  section  2.2  have  the  disadvantage  that 
they  give  energy  release  rates  or  related  quantities  but  explicit 
determination  of  stress  intensity  factors  may  not  be  possible  when  mixed 
mode  situations  are  encountered.  It  should  be  noted  however  that  the  general 
approach  of  section  (2.2)  (equations  (8.33)  to  (8.36))  is  not  restricted  to 
linear  stress-strain  laws.  An  alternative  approach  has  been  developed  by 
Stern  and  co-workers  (20]  to  [24]  who  derive  invariant  integrals  for  plane 
linear  elastic  problems  by  means  of  Betti's  reciprocal  theorem. 

The  starting  point  of  their  analysis  is  the  reciprocal  theorem  written  as 


/  “  T.u)  ds  =  0 


(8.34) 


where  3R  is  the  boundary  of  a  plane  simply  connected  bounded  region  R  . 

The  states  (a...  u.)  and  (a..,  u.)  are  two  distinct  equilibrium  states 
ij  i  ij  i 

(o„,  °i^  ^eing  the  state  corresponding  to  a  given  boundary  value  problem 
the  other  state  is  an  auxiliary  one;  T  and  T  are  the  boundary  tractions 
associated  with  these  elastic  states.  For  plane  crack  problems  a  contour 

3R  such  as  that  shown  in  Fig.  4  is  taken. 

{  i  .  <•  u 

The  idea  of  the  method  is  to  choose  the  auxiliary  state  so  that  the 

integral  around  C  gives  the  coefficient  of  the  required  stress-singular i tv 
e 

as  e  tends  to  zero.  For  a  crack  problem  if  the  auxiliary  solution  also 
satisfies  the  6tress  free  boundary  condition  along  the  crack  it  is  possible 
to  relate  the  integral  around  C£  to  the  integral  around  by  using 
equation  (8.54).  The  result  is  that 
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lim  J  (  T .  u  -  T.u)ds 
e+0  C 


J  (T.fi  -  T . u )  ds 
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(8.55) 


Now  the  stresses  and  displacements  in  the  neighbourhood  of  the  crack 

tip  referred  to  the  natural  polar  co-ordinate  system  shown  in  Fig.  4  are:- 

»r  lh VafD»-OC«£  -C.*»S3lC(  -L^X +s(r"V 
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where  u®  and  u^  are  the  radial  and  tangential  components  of  the  displacement 
u°  of  the  crack  tip  and 
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(8.57) 


are  the  conventional  stress  intensity  factors.  The  remainder  terms  are  of 


the  order  indicated  in  distance  from  the  crack  tip. 

Thus  when  C  is  a  small  circular  contour  close  to  the  crack  tip, 
ds  “  rd8,  and  the  dominant  contribution  from  the  traction  T  involves  tin 


stress  intensity  factors  multiplying  spatially  varying  terms  proportional 
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to  r 


*/2 


Hence  if  the  auxiliary  displacement  £  can  be  chosen  to 


-l/: 


have  dominant  behavior  proportional  to  r  as  r  tends  to  zero  the  product 


of  T.u  with  ds  will  have  a  non-zero  finite  limit  as  e  tends  to  zero. 


3/2 


Similarly  the  traction  T  which  will  now  be  proportional  to  r  as  r 


tends  to  zero  will  lead  to  the  required  1/r  behavior  in  the  product 
0 


T.(u-u  ).  The  only  stipulation  is  that  u  and  T  must  be  equilibrium 
solutions  and  moreover  must  satisfy  the  stress  free  crack  boundary 
conditions.  For  elastic  problems  an  exact  solution  of  the  field  equations 
satisfying  stress  free  boundary  conditions  can  be  found  by  an  eigenfunction 
approach.  The  auxiliary  elastic  field  required  for  the  above  problem  has 
been  given  by  Stern  et  al  [20]  as:- 


u  =  - r -  {[(2k  +  1)cos  -  3  cos  ^  )c  +  [(2ic+l)sin  mr-  -  sin  -  ]c  1 

2(2*r)*(l+«)  1  222 


u  =  - r - {  [-( 2<-l  )s  in  ^  +  3  sin  -|]  c  +  [  ( 2*r- 1  )cos-^  -  cos  ^]c_) 
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- -  {  [7  cos  ^  -  3  cos  f  ]  c  +  [7  sin  ^  -  sin  f  ]c,} 
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- - - ([cos  +  3  cos  |-]c  +  [sin  ^  +  sin  -lit.1  (8.38) 
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or  o  = - - { 3  [sin  ^  +  sin  |]c.  -  [cos  ^y-  +  cos  |  Jc  } 
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where  c^  and  C2  are  arbitrary  constants. 

Now  on  the  inner  circular  boundary  the  evaluation  of  the  contour 


integral  in  terms  of  the  traction  and  displacement  (relative  to  that  of 
the  crack  tip)  components  takes  the  form 


I  *  -  f  [(u  -  u  )-t  -  u-t]ds 

c. . 
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Upon  substitution  from  (8.56)  and  (8.58),  a  routine  evaluation  of  the 
integral  produces 
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(8.59) 


where  the  remainder  term  goes  to  zero  with  t  as  indicated.  Thus,  for 
arbitrarily  small  (  (8.55)  produces  the  representation  formula 


:iKi +  c2kh  =  icI(-'-0)'-  ■  i’y* 


(8.60) 


where  it  is  important  to  note  that  the  contour  C  involves  only  the  outer 

boundary  since  both  t  and  ^necessarily  vanish  on  the  crack  faces. 

Furthermore,  the  rigid  body  displacement  may  be  discarded  in  the 

0  - 

evaluation  of  (8.60)  since  the  contribution  of  u  't  on  any  closed  contour 
vanishes  even  when  the  origin  is  contained  in  the  interior.  It  remains  only 
to  obtain  u  and  t_  on  the  outer  boundary  from  the  prescribed  data  so  that 
the  contour  integral  may  be  evaluated  as  a  linear  combination  of  c^  and  Cj, 
the  coef ficienUof  which  are  the  desired  stress  intensity  factors  and  K  • 
The  idea  of  the  above  method  can  in  principle  be  extended  to  a  number 
of  other  situations,  however  a  major  restriction  is  that  the  auxiliary 
solution  must  be  a  full  solution  of  the  equilibrium  equations.  Nevertheless, 
for  a  variety  of  elastic  problems  sue  h solutions  can  be  found  and  Stern  and 
co-workers  have  applied  the  method  to  a  number  of  other  problems  f ( 2 1 1  to  [2^]) 
An  obvious  advantage  of  the  method  is  of  course  the  computation  of  near  iield 
crack  tip  information  from  an  integral  evaluated  far  from  the  crack  tip. 


A  variety  of  integral  equation  methods  have  been  suggested  in  order 
to  overcome  the  difficulties  outlined  in  the  introduction.  However,  most 
of  these  methods  have  been  designed  to  treat  a  special  class  of  problems 
and  are  not  easily  generalised.  In  the  following  subsections  these 
methods  are  outlined. 

We  begin  by  briefly  discussing  how  the  integral  equation  (8.15)  can 
be  used  for  the  fracture  problem.  In  section  3.1  it  is  shown  how  symmetry 
can  be  used,  in  certain  cases,  to  reduce  the  problem  to  a  boundary  value 
problem  for  which  the  formulation  discussed  in  the  introduction  (equation 

(8.15) )  is  appropriate.  In  section  3.2  partitioning,  dividing  the  body 
into  regions,  is  used  together  with  continuity  conditions  along  an 
extension  of  the  crack  boundary  to  couple  the  two  regions  together.  In 
both  these  methods  integral  equations  like  equation  (8.15)  are  used  so  no 
special  attention  is  paid  to  the  crack  as  far  as  the  integral  equation 
formulation  is  concerned.  In  section  4  we  discuss  various  ways  of 
approximating  the  integral  equation  in  order  to  accurately  model  the  crack. 

In  sections  3. 3-3.6  we  discuss  integral  equation  formulations  which 
pay  particular  attention  to  the  crack  geometry.  In  section  3.3  a  method 
which  uses  exact  'source-crack'  Green’s  functions  is  outlined.  The 
limitation  of  this  method  is  that,  except  for  special  crack  configurations 
(e.g.  straight  cracks  in  plane  strain  or  stress),  very  few  of  these  Green's 
functions  are  known.  In  section  3.4  straight  cracks  modelled  by  Volterra 
dislocations  are  considered,  thus  an  integral  equation  is  obtained  for 
the  crack  in  terms  of  a  continuous  distribution  of  (virtual)  dislocati  mis 
This  is  coupled  to  a  boundary  integral  formulation  such  as  that  of  equation 

(8.15)  in  order  to  deal  with  finite  specimens.  The  limitation  of  this 
approach  is  that  only  plane  problems  and  straight  cracks  are  envisaged. 
Nevertheless,  it  is  natural  to  represent  the  crack  by  elements  which 


represent  a  jump  in  displacement  and  the  dislocation  is  a  fundamental 
unit  of  this  kind.  More  fundamental  than  a  Volterra  dislocation,  however, 
is  the  Somigliana  dislocation  which  is  an  infinitesimal  dislocation  of 
area  dS  with  an  associated  discontinuity  vector  (representing  the  jump 
in  displacement  Uj  across  dS)  which  is  effectively  constant  over  the  small 
area  dS.  We  discuss  crack  modelling  by  Somigliana  dislocations  in  Section 
3.5.  Finally,  in  section  3.6  we  show  how  the  dislocation  models  of 
sections  3.4  and  3.5  can  be  derived  naturally  from  Betti's  reciprocal 
theorem. 


3 . 1  Use  of  symmetry  (where  possible)  to  replace  crack  by  boundary  value 
problem 

<  W  " 

For  certain  symmetric  geometries  and  loading  conditions  it  is  possible 
to  replace  a  crack  problem  by  a  mixed  boundary  value  problem  and  thus 
apply  the  standard  boundary  integral  equation  of  equation  (8.15).  We 
illustrate  this  approach  by  considering  a  fracture  specimen  such  as  shown 
in  Fig.  5.  If  the  applied  loading  is  tensile  and  symmetric  with  respect 
to  Xj  ■  0,  it  is  clear  that  fracture  should  occur  in  the  opening  mode  and 
the  stress  analysis  of  the  configuration  (Fig.  5)  will  lead  to  a  mode  I 
stress  intensity  factor.  Using  the  symmetry  of  the  specimen  it  is  sufficient 


to  analyse  the 

region 
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>  0. 

For  a 

traction 

free  crack  the  boundary 

conditions  on 
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on 

x2 

■=  0 

Xj  6  AD 

t2  -  0 

on 

x2 

=  0 

xL  €  AO 

(8 

u2  ■  0 

on 

x2 

=  0 

x 1  6  OD  . 

The  usual  notation  for  stress  and  displacement  has  been  used,  e.g.  (uj.u^) 
are  the  displacements  in  the  (x(,x2)  directions,  and  the  traction  t.  =  ’..n., 
n-  being  the  components  of  the  outward  drawn  normal.  We  write  the  integral 


equation  (8.15)  in  the  equivalent  form 


/  [u. (Q)  -  u.(P)]T..(P,Q)dS(Q) 

s  1  1  J1 

=  J  t . (Q)lK . (P,Q)dS(Q) 


(8.62) 


where  S  is  the  boundary  ABCD ,  i  =  1,2,  and  the  integral  equation  is 

to  be  satisfied  for  all  P  €  S.  With  either  the  traction  vector  t.  or 

i 

displacement  vector  u^  specified  on  the  boundary  ABCD,  consistent  with 
the  assumed  symmetry,  the  above  integral  equation  can  be  reduced  to  a 
system  of  simultaneous  equations  as  discussed  in  section  A  or  elsewhere 
in  this  series.  Of  course  accurate  element  modelling  is  required  to  obtain 
good  results  for  the  singular  crack  tip  stresses  etc.  Note  that  the  integral 
equation  is  required  along  the  whole  of  boundary  AD  whereas  the  attempt  to 
directly  model  the  crack  boundary  discussed  in  the  introduction  would 
have  led  (if  successful)  to  an  integral  equation  only  along  the  crack  AO 
in  addition  to  the  external  boundary. 

For  some  early  work  on  applications  of  this  method  to  fracture 
specimens  where  a  symmetry  plane  allows  the  problem  to  be  set  up  as  a 
boundary  value  problem  see,  e.g.  Cruse  and  Van  Buren  [25].  In  [25]  the 
thickness  of  the  fracture  specimen  was  taken  into  account,  the  boundary 
and  crack  surfaces  being  divided  into  triangles  in  which  the  displacements 
and  tractions  are  assumed  constant.  In  section  A  we  show,  for  an  anti¬ 
plane  elastic  problem,  how  more  refined  element  modelling  can  be  effective. 


3 . 2  Partitioning  (subdividing  the  body  into  distinc t  regie ns) 

I 

When  there  is  no  symmetry  in  the  cracked  specimen  it  is  still 
possible  to  get  a  formulation  like  that  obtained  in  the  last  section  bv 
subdividing  the  body  into  regions.  We  illustrate  this  approach  bv  considering 
the  configuration  shown  in  Figure  6.  For  convenience  ••e  denote  the  part 


of  the  boundary  described  by  ABCDA  as  and  the  part  described  bv  ADF.FA  as  V  ^ 


Then  applying  the  integral  equation  (8.62)  to  both  regions  (i)  and  (ii) 
gives  the  equations 


J  Iu.(1)(Q)  -  u.(1)(P)]T. .(P.Q)dS(Q) 
r  1  1  J1 

1 

=  J  t.(1)(Q)U..(P,Q)dS(Q) 

r  1  J1 

1 

and 

/  [u.(2)(Q)  -  u.(2)(P)]T..(P,Q)dS(Q) 

r2 

=  J  t.(2)(Q)U..(P,Q)dS(Q) 

r  1  Ji 


(8.63) 


(8.64) 


Along  the  common  boundary  AD  we  have  the  boundary  conditions 


(i)  Continuity  of  stress,  i.e.  t.^^ 

(ii)  Traction  free  crack  condition  t. 

i 

(iii)  Continuity  of  displacement  on  OD, 


(2)  c 

=  -t^  for  all  points 

=  0  on  AO. 

(1)  (2) 
u  .  =  u  .  on  05. 

l  l 


on  AD . 


It  is  possible  to  treat  the  above  equations  in  the  usual  way  although, 
of  course,  accurate  modelling  is  required  in  order  to  obtain  reliable 
results.  For  an  early  application  of  this  approach  see,  e.g.  Rizzo  and 
Shippy  [26]  where  a  formulation  of  this  kind  is  given  for  non-homogeneous 
(two  dimensional)  inclusion  problems. 


3 . 3  The  use  of  exact  'source-crack'  Green's  functions 

One  way  of  avoiding  the  difficulty  discussed  in  the  introduction  is 
to  avoid  the  crack  altogether  by  using  as  a  source  function  in  the  integral 
equation  formulation  the  solution  to  the  problem  of  a  point  source,  placed 
anywhere,  and  a  crack  in  an  infinite  medium.  This  approach  thus  relies  on 
the  availability  of  such  a  solution  and  requires  that  the  solution  be 


i 


1 


known  at  all  points  of  the  body,  i.e.  in  particular  at  the  boundary  of 
the  body  where  the  integral  equation  is  applied.  These  requirements 
restrict  the  applicability  of  such  an  approach  since  solutions  are  known 
only  for  special  crack  geometries.  In  three  d imensions _for  example  a 
crack  with  a  flat  elliptical  planform  is  about  the  limit  of  cases  that  can 
be  treated  by  exact  analysis  and  the  stress  field  everywhere  takes  on  a 
very  complicated  form.  Some  applications  of  this  method  to  problems  of 
plane  anisotropic  elasticity  have  been  made  by  Snyder  and  Cruse  [27],  The 
resulting  integral  equation  on  the  boundary  B  takes  the  form 


u.(z  )/2  +  j  T'.(z  ,z  )u.(z  )dS( z  )  =  /  U'.Cz  ,z  )t.(z  )dS(z  ) 
JO  jikOik  k  J  ji  k  0  l  k  k 


(8.65) 


where  z  =  x  +  ix  is  a  point  on  the  boundary  B ,  z,  =  x,,  +  ix„,  ,  and 

0  10  20  k  lk  2k 

the  integrals  are  to  be  interpreted  in  the  principal  value  sense.  The 

displacements  and  tractions  U .  and  T'. .  are  derived  from  the  fundamental 

J1  J1 

solution  for  the  infinite  plate  with  a  crack  lying  on  Xj  =  0  -a  <  <  a 

and  with  unit  loads  applied  in  the  x^  direction  acting  at  some  point 
^X10’X20^  not  !yin8  on  t*le  crack.  The  functions  Ul^  and  T^  given  in  [27] 
are  for  orthotropic  materials.  Note  however  the  corrections  indicated  in 
the  later  paper  by  Cruse  [5] .  Thus  in  their  applications  only  mode  I  and  II 
stress  intensity  factors  are  found.  It  should  be  noted  however  that  for 
general  anisotropic  solids,  plane  elastic  problems  can  result  in  stress 
intensity  factors  of  modes  I,  II  and  III  even  though  the  displacement  field 
does  not  vary  in  the  three  direction.  The  method  of  this  section  can  be 
extended  to  this  situation  also,  see  the  papers  [28,29,30]  for  information 
on  the  necessary  fundamental  solutions. 

As  discussed  earlier  it  is  only  for  certain  special  problems  that  the 
fundamental  solution  of  the  unit  load-crack  interaction  is  known  in  a 
suitably  compact  form  and  thus  this  method  is  of  limited  application. 


I 


Nevertheless,  for  completeness  we  briefly  sketch  the  derivation  of 
equation  (8.65)  following  [27],  Betti's  reciprocal  theorem  in  the  absence 
of  body  forces  can  be  written 


/  T ! . u . d  S  =  f  Ul.t.dS  (8.66) 

B+L+r  J1  1  B+L+T  J1  1 


where  L  is  the  crack  boundary,  B  the  boundary  of  the  body,  and  7  is  a 
circle  of  vanishing  radius  e  surrounding  the  unit  load  applied  at  the 
point  (x^.x^q).  Since  the  crack  L  is  assumed  to  be  traction  free  the 
integrals  in  equation  (8.66)  over  the  crack  L  are  both  identically  zero,  i.e. 


/ 

L 


T! .u.dS 
Jl  1 


/  u:.t.ds  =  o 

L  1 


(8.67) 


The  boundary  integral  equation  (8.65)  is  thus  derived  in  the  usual  manner 
by  letting  V  ■*  0  and  (xj0,X20^ 

In  [31]  a  somewhat  similar  situation  is  considered. 


3 .4  Integral  equations  for  ideal  crack  geometries 

In  the  last  section  the  crack  problem  was  treated  by  using,  as  a 
fundamental  solution,  the  complete  stress  and  displacement  field  due  to  a 
unit  force  interacting  with  a  crack.  However,  as  indicated  in  that 
section  such  solutions  are  known  only  for  a  special  class  of  crack  geometries 
e.g.  straight  cracks  in  plane  strain  configurations.  The  next  stage  in 
versatility  is  perhaps  reached  by  modelling  the  crack  by  extended 
dislocations  (Volterra  dislocations).  Crack  problems  formulated  in  this 
way  have  been  reviewed  by  Bilby  and  F.shelby  [32]  for  the  homogeneous  medium 
case,  and  integral  equations  for  cracks  in  bimaterials  have  been  derived 
and  treated  numerically  in  Atkinson  [33]  and  Cook  and  Erdogan  [ 3 '• )  .  This 
approach  has  been  coupled,  with  the  boundary  integral  equation  method  (8.15), 
by  Rudolph  and  Ashbaugh  [35].  They  treat  the  two  dimensional  problem  of  a 


» 


finite  length  crack  in  a  bounded  linearly  elastic  isotropic  medium 
subjected  to  in  plane  forces.  The  superposition  principle  is  used  to 
couple  together  the  B.l.E.  method  applied  to  the  unflawed  medium  with 
boundary  stresses  those  induced  on  the  boundary  by  the  'perturbed  problem', 
that  of  a  crack  in  an  infinite  medium.  Since  this  'perturbed  problem'  must 
represent  a  stress  free  crack,  the  traction  on  the  crack  surface  must  be 
zero  thus  producing  an  integral  equation  along  the  crack  surface  which 
involves  the  interior  stresses  induced  at  the  crack  boundary  by  the  B.l.E. 
solution,  hence  coupling  the  two  problems. 

The  treatment  given  in  [35]  concentrates  on  a  mode  1  problem  (hence 
a  certain  amount  of  symmetry  is  assumed),  although  in  principle  the  method 
can  be  applied  (to  straight  cracks)  when  there  are  other  modes  present 
(such  a  situation  is  considered  by  Chang  and  Morgan  [36]).  The  stress 
field  of  a  mode  I  crack  in  an  infinite  medium  is  expressed  in  terms  of  a 
density  of  virtual  edge  dislocations  f  (  O  in  the  form 


(8.68) 


'  i^TT.{  fC«K..<s.OaE 


where 


Kj  =  (k-1)  In  r^  -  Ax^/r^ 

C-x.  Ax2(E-x1) 

K,  «  -2(^+1  )arctan  ( - -)  -  - ^ - 

2  x2  r2 


(8.69) 


2  2  2 

r  *  (5-*, )  ♦  x„  and  the  kernels  K. .  can  be  derived  from  the  K.  by 

12  ij  i 

differentiation.  Note  that  the  dislocation  density  can  be  expressed  in 


terms  of  the  u2  displacement  by 


JO 


,,  v  3  r  (2),  +*_  3  ,  (2). 

f(xl)  37^  rU2  (xl'°  )]  =  "  3^  U2  (xl'°  )] 


(8.70) 


where  the  crack  is  taken  to  lie  on  Xj  =  0,  -1  <  x^  <  1 . 

The  solutions  for  the  'unflawed  medium'  problem  (designed  by  a 
superscript  (1)  in  our  notation)  and  the  crack  in  an  infinite  medium 
(designated  by  superscript  (2))  are  now  combined  to  obtain  the  solution 


for  the  original  problem.  Let  denote  that  portion  of  C  where  u^  is 


specified  and  C. 


'O 


^denote  that  portion  where  t^  is  specified,  then 


c  -  c.u  ♦  c.c 

1  1 


i  =  1,2 


and 


(8.71) 


u.(1)(P)  +  u/2)(P)  *  u.(P),  P  e  c.u 

ill  l 


t.(1)(P)  ♦  c.(2)(P)  -  c.(p>’  p6Ci 

11  1 


where  u^(P)  and  t^(P)  are  the  specified  boundary  displacements  and 
tractions  of  the  original  problem. 

Also  for  |x^|  <1,  x2  *  0 

(1)  (2)  (8.72) 

t2  tx^.O)  +  t2  ;(Xj,0)  =  0 

t1(1)(x1>0)  =  0 

the  last  condition  is  an  im^plied  restriction  on  the  symmetry  of  the 
boundary  C  and  loading  system  such  that  only  one  fracture  mode  is  operative 
(mode  I).  [Note  of  course  that 


t2(x1>0)  =  o22(xj,0)  and  t^Cx^.O)  =  o.^x^.O)  ) 

The  use  of  the  above  conditions  leads  to  a  pair  of  coupled  integral 
equations,  the  integral  over  C  being  of  the  usual  B.I.E.  type.  The  integral 
equation  along  the  crack  has  a  Cauchy  integral  part  due  to  the  dislocation 
representation  of  the  crack.  Thus  from  equation  (8.68)  the  normal  stress 


on  the  crack  line  reduces  to 


.  (2),  _  (2),  ^  4U  f  f(c)d( 

2  (xl’0)  "  °22  (xl*0)  “  „(<+!)  c-xj 


If  the  crack  is  to  be  closed  at  both  ends,  the  subsidiary  conditic 


(8.73) 


J  f(Od£  =.  0 


must  be  satisfied.  A  well  known  inversion  formula  for  the  integral 
equation  (8.73)  then  gives 

Ait(l-x  2)%f(*  )  =  -  /  fp-Q-  t  (2)U,0)d£  +  d 

11  ii  (C-x^  2 

for  | x ^ |  <1 


(8.74) 


(8.75) 


where  d  is  an  arbitrary  constant  to  be  determined  from  condition  (8.74). 
Note  that  the  above  singular  integrals  are  to  be  interpreted  as  their 
Cauchy  principal  values. 

Since  from  the  boundary  condition  (8.72) 


t2(2)(£,0)  =  -t2(1)(C,0)  , 


kl  <  i 


(8.76) 


the  integral  equation  (8.75)  contains  double  integrals  (t2^(£,0)  =  V,n) 
i6  given  for  the  unflawed  medium  by  an  equation  like  equation  (8.14)  of 
the  introduction).  Fudolph  and  Ashbaugh  [35]  overcome  this  difficulty  by 
changing  the  order  of  integration  and  evaluating  various  singular  integrals 
analytical ly . 

Finally  the  integral  equations  (8.13)  and  (8.75)  are  solved  numerically 

for  the  unknown  density  f(x,),  -1  <  x,  <  1  and  the  unknown  u.(P)  and  t . < F ) 

1  1  J  1 

on  corresponding  parts  of  the  boundary  C.  It  should  be  noted  that  the 
efficacy  of  the  procedure  depends  to  a  certain  extent  on  quite  a  bit  of 
subsidiary  analytical  work  in  reducing  various  double  integrals  to  single 


•.  -.  „\ 


3.’ 


ones.  An  alternative  to  the  step  from  (8.73)  to  (8.75)  would  be  to 
solve  numerically  the  integral  equation  (8.73)  directly  with  appropriate 
modelling  of  f(£)  • 


3.5  Crack  modelling  by  Somigliana  dislocations 


!  \  •  ^  t  <*  ... 

A  natural  extension  to  the  approach  of  the  last  section  is  to  model 


the  crack  as  a  mosaic  of  infinitesimal  dislocations  of  area  dS„ ,  normal  n, 

1  k 

and  a  discontinuity  vector  which  is  effectively  constant  all  over  the 
small  area  dS^.  The  displacement  produced  by  one  of  these  elementary 
dislocations  can  be  written 


du.(r)  -  b.n  dS.p..^  (r  -  r') 
J-  lklik- 


(8.77) 


where  p^  (r  -  r')  is  the  stress  produced  at  a  point  r’  on  by  a  point 
force  of  unit  magnitude  at  r  parallel  to  the  x^  axis  and  n^  is  the  normal 
to  at  r1.  A  Somigliana  dislocation  at  (see  Fig.  7)  is  characterised 

by  the  vector  b  which  is  the  displacement  on  the  arrow  side  of  defined 
by  njminus  the  value  on  the  tail  side  of  the  arrow.  Thus  in  an  infinite 
linear  elastic  medium  the  displacement  field  due  to  the  Somigliana  dislocation 
can  be  shown  to  be 


Uj(')  “  f  bi(-')pik)(-  ”  nkdSl 


(8.78) 


Note  that  as  far  as  the  above  formulation  is  concerned  the  surface  S^ 
can  take  any  shape,  so  in  principle  a  crack  of  any  shape  can  be  represented 
in  this  way.  However,  the  numerical  solution  of  the  integral  equation  which 
results  when  expression  (8.78)  is  differentiated  to  obtain  the  stress  field 
requires  further  investigation.  A  proof  of  (8.78)  which  also  verified  its 
properties  for  an  anisotropic  medium  was  given  by  Burgers  [37).  A  derivation 
of  (8.78)  by  means  of  the  reciprocal  theorem  is  attributed  to  Esholby  in 


3.1 


the  book  by  Nabarro  [38].  We  shall  give  in  the  next  section  a  derivation 
based  on  Betti's  reciprocal  theorem. 

In  principle  the  integral  equation  deduced  from  (8.78)  by  computing 
the  stresses  at  r  and  then  letting  r  tend  to  can  be  coupled  with  a 
boundary  integral  formulation  of  the  external  boundary  in  a  similar 
manner  to  that  of  the  last  section.  However,  the  numerical  implementation 
of  such  a  procedure  has  yet  to  be  made.  It  should  be  noted  nevertheless 
that  this  formulation  can  be  used  for  any  crack  profile. 

The  special  case  of  a  flat  crack  in  an  infinite  medium  has  been 
tackled  by  Weaver  [39]  using  a  formulation  similar  to  that  discussed 
above.  He  reduces  the  singular  nature  of  the  integral  equations  by 
first  doing  an  integration  by  parts. 

3 . 6  General  formulation 

We  discuss  here  a  general  formulation  of  the  crack  problem  in  which 
the  crack  profile  is  modelled  by  a  distribution  of  Somigliana  dislocations. 
We  begin  with  the  formulation  of  the  situation  shown  in  Fig.  1  of  section  1. 
If  in  equation  (8.19)  the  external  boundary  S  is  allowed  to  tend  to  infinity 
and  the  boundary  displacement  and  stresses  are  such  that  the  integrals 
over  S  tend  to  zero,  then  since  It^(Q)  =  0  on  r  the  following  result  is 
obtained 


It  is  easily  seen  that  this  result  is  equivalent  to  that  of  equation  (8.78) 

by  noting  that  T..  *  -p.^^n  if  n  is  defined  as  in  Fig.  7.  Note  that 

in  the  formulation  of  Weaver  C.39^  where  a  crack  with  a  flat  three  dimensional 

plan  form  is  considered,  the  surface  r  is  associated  with  the  upper  crack 

+ 

surface  and  hence  has  unit  normal  n.  *  -6 . _  in  his  notation,  i.o.  it  is  in 


the  opposite  direction  to  that  of  Fig.  7.  As  indicated  in  the  last 
section,  the  stresses  can  be  computed  from  (8.79)  by  differentiating  with 
respect  to  the  point  p. 


In  general  for  an  internal  crack  such  as  shown  in  Fig .  1  we  have 
(equation  (8.19)) 

(Q)dS(Q)  +  J  U.  .(p,Q)t.(Q)dS(Q) 

S  1J  J 


u.(p)  =  -  J  T..(p,Q)u. 

S  J  J 


-  /  T. j(p,Q)AUj(Q)dS(Q) 


(8.80) 


where  Au.(Q)  =  u.(Q  )  -  u.(Q  )  and  we  have  used  the  fact  that 
J  J  J 

Itj(Q)  =  t ^(Q  )  +  t j ( Q~ )  =  0  across  the  crack.  Differentiating  (8.80)  with 
respect  to  the  co-ordinates  of  p  to  obtain  the  stresses  gives 


a . . (p) 

iJ 


(p,Q)tk(Q)dS(Q)  -  / 
S 


Skii(p’Q)uk(Q)dS(Q) 


-  J  S  .(p,Q)Au  (Q)dS(Q) 

r  k!  J  K 


(8.81) 


The  functions  D,  .  .  and  S,  .  .  are  obtained  by  differentiating  U.  .  and  T. . 

kij  kij  *  ij  U 

respectively,  the  latter  two  functions  can  be  deduced  from  equations  (8.8) 
and  (8.7)  of  the  introduction. 

There  is  no  difficulty  with  the  integral  equation  (8.80)  as  p  tends 

to  a  boundary  point  P  of  S.  The  result  is  an  integral  equation  like  (8.80) 

with  p  replaced  by  P  and  u^(p)  on  the  left  hand  side  of  (8.80)  replaced  by 
u^(P)/2  (cf.  equation  (8.15)  of  the  introduction).  Thus  for  P  belonging 
to  S  the  usual  boundary  integral  equation  applies  with  the  addition  of  a 
non-singular  integral  involving  the  unknown  Au^(Q).  When  p  tends  to  a 
point  of  T,  we  have  already  seen  in  the  introduction  that  the  formulation 
(8.80)  is  inadequate.  However,  in  our  opinion,  there  is  nothing  wrong  with 
an  integral  equation  based  on  (8.81)  for  points  p  tending  to  P  belonging  to  ! 


The  traction  free  (or  stress  free)  houndary  condition  for  a  crack  is 
sufficient  to  specify  the  integral  equation  for  Au^(Q),  Q  P  V. 

The  difficulty,  of  course,  is  in  the  implementation  of  this  formulation. 
Nevertheless,  we  feel  it  shoulJ  even  be  possible  to  deal  numerically  with 
the  improper  integrals  resulting  in  (8.81)  when  p  -*■  P  6  f.  Failing  this, 
of  course,  it  should  be  possible  to  integrate  by  parts  in  the  integral 
over  F  in  (8.81)  before  taking  the  limit  p  -*•  P  6  r .  This  was  done  in  [34] 
for  the  case  of  an  infinite  medium.  Such  a  procedure  does  of  course 
reduce  the  singular  nature  of  the  integrals  and  for  two  dimensional  problems 
it  leads  to  singular  (Cauchy)  integrals  such  as  considered  in  section  3.4. 
This  procedure  is  analogous  to  modelling  the  crack  by  extended  dislocations 
(Volterra  ones)  as  opposed  to  the  infinitesimal  Somigliana  dislocations. 

An  integral  equation  formulation  somewhat  similar  to  that  of  (8.81) 
for  the  crack  has  been  formulated  and  implemented  by  Montulli  [ 40 )  . 

However,  he  uses  such  an  integral  equation  together  with  a  standard 
boundary  integral  equation  along  the  crack  surfaces.  A  somewhat  similar 
procedure  is  being  developed  by  Xanthis  [41].  In  either  formulation  we 
do  not  see  that  the  extra  equations  are  necessary,  all  that  is  required  to 
model  the  crack,  in  our  opinion,  is  the  dislocation  distribution.  It  should 
be  stressed  that  the  above  formulation  is  intended  as  a  general  one  to 
deal  with  any  crack  shape,  the  success  or  otherwise  of  this  approach  depends 
of  course  on  the  implementation.  There  are  other  methods  which  have  been 
used  in  two  dimensions  to  deal  with  notches  and  corners;  we  mention 
in  this  respect  the  work  of  Barone  and  Robinson  [42] .  Clearly  more  work  is 
required  before  the  definitive  approach  is  determined. 


4 .  Modelling  and  numerical  results 

We  have  indicated  throughout  the  earlier  sections,  and  in  the 
introduction,  that  precise  modelling  of  the  crack  tip  singularity  is 
essential  to  an  accurate  analysis  of  the  fracture  problem.  In  this 
section  we  give  an  account  of  some  recent  work  (Xanthis  et  al  [43 j) 
which  pays  due  attention  to  this  singularity.  The  situation  considered, 
that  of  anti-plane  strain  deformation  (longitudinal  shear),  is  perhaps  the 
simplest  in  which  this  singular  behavior  can  be  exhibited.  Nevertheless, 
the  ideas  discussed  here  should  be  of  use  in  the  more  complicated  problems 
discussed  in  section  3. 


4.1  General  formulation 


Since  anti-plane  strain  deformation  is  being  considered  there  is 
only  one  non-zero  displacement  u  and  this  acts  perpendicular  to  the  plane 


of  the  paper.  The  general  problem  is  to  determine  u  in  a  bounded  open 
2 


region  12  in  R  with  boundary  312  where  u  satisfies  the  equations  and  boundarv 
conditions: 


2 

V  u  =  0 

in 

12 

u  ■  u 

on 

312, 

3u/3h  =  t 

on 

312, 

3u/3n  +  cu  =  d 

on 

312 

(8.82) 


312 


j,  3122  and  3 Q ^  are  disjoint  portions  of  the  boundary  with  union  equal 


to  312  i.e.  312  =*  312 1  U  3122  \J  3123  and  312i  f\  312  ^  =  <f  i  t  j .  The 
expressions  u,  t,  c  and  d  are  given  functions  and  3u/3n  is  the  outward 
normal  derivative  of  u).  Since  the  displacement  u  satisfies  Laplace's 
equation  a  scalar  version  of  equation  (8.15)  can  be  derived  via  Green's 
third  identity.  The  resulting  integral  equation  on  the  boundary  31  (  ol  the 
body  is 


j  T(P,Q)lu(Q)  -  u(P)]dS(Q)  =  j  U(I’.Q)  (Q)dS(Q) 


(8.83) 


where 


U(P,Q)  3  (  l/2n  )ln( 1/r) 


for  all  P  €  Si;  , 


with  r  =  | P-Q I 


T(P,Q)  =  VqU(P ,Q) 'n(Q)  where  n.(Q)  is  the  outward  unit  normal  at  the 
point  Q  €  SSI  and 

=  7u(Q)  .  n(Q)  . 

dn  — 


For  Q  6  S.Q  let 


u*(Q)  3  u(Q),  t  * (Q)  3  P-  (Q) 

3n 


then  u*  and  t*  satisfy 


f T ( P , Q ) [u*(Q)  -  u* (P) ] dS(Q )  =  Ju(P,Q)t*(Q)dS(Q) 


u*(P)  =  u(P) 

t*(P)  =  t(P) 


V  P  6  SP 


v  p  e  Si 


v  p  e  jn. 


(8.84) 


(8.85) 


t*(P)  +  c(P)u*(P)  =  d(P)  V  P  6  Si), 


The  object  is  thus  to  compute  approximations  u  and  t  to  u*  and  t*  from 
the  integral  equation  (8.83)  and  boundary  conditions  (8.85). 

Any  particular  boundary  element  method  is  characterised  by: 

(a)  the  choice  of  the  classes  of  functions  used  to  approximate  u*  and  t*, 
and  (b)  the  way  in  which  the  integral  equation  (8.83)  is  used  in  the 
determination  of  u  and  t.  In  Xanthis  et  al  [43]  u*  and  t*  are  approximated 
by  linear  combinations  of  computationally  convenient  functions  \\U\ 


=  1 . n1  ,  and  v  ,  ,  £'  =  1 . n~  ,  respectively,  i.e.  one  writer 
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it 


(e  3) 

and  P  ’  .  Tne  approximat ion  tc  the  part  of  the  boundary  passing 

(e) 

through  these  nodes  is  specified  by  P  (O,  (,  6  10,1]  where 


p(e)(o 


i  =  l 


N.  (2)(C)P(e  ,:l) 
i 


(£.69) 


with 


V 


N  (2)U)  =  (1-0(1-20.  k2(2)(c)  =901-0,  k3(2)(0  =  0201) 

(2) 

The  K.  ,  i  =  1,2,3  are  the  quadratic  Lagrange  polynomials  defined  on 


(2)  (2) 

[0,11  with  the  ore  pert  y  N .  (O  )  =  6..  where  Z, '  =  (  i-1 ) /2  , 

i  j  ij  j  - 


(2) 


j  =  1,2,3. 


(p) 


More  generally  one  can  define  K.  ,  i  =  l,...,p+l  to  be  the  p 


th 


order  Lagrange  polynomials  with 


s.'Oe.'O  •  c.lp)  ■  (i-D/p, 

i  '1  1 J  J 


j  =  1 , • • • ,P+1 


For  a  good  approximation  to  the  boundary  the  element  boundary  nodes  need 
to  be  chosen  so  that  the  portion  of  the  boundary  between  them  is  smooth, 
hence  corners  of  of.  are  chosen  as  element  boundarv  nodes. 


A . 1 . 2  Choice  of  approximating  functions 

The  functions  v^U"  and  arc  chosen  with  the  following  considerations 

in  mind:  (i)  u*  is  continuous  but  t*  may  have  discontinuities,  e.g.  at 
corners  or  boundary  points  where  there  is  a  change  in  boundary  condition 
type.  Thus  the  v^U^'s  are  chosen  to  be  continuous  and  the  v^^'s  such  that 
functions  with  discontinuities  at  selected  boundary  points  can  be  represented 
accurately,  (ii)  the  functions  v^u  and  must  be  computational Iv 

convenient . 

It  is  assumed  that  the  element  nodes  have  been  chosen  so  that  the 


mapping  (£.89)  is  invertible  so  an  inverse  mappin; 


£  =  f(c)(P(£)(£))  , 


can  be  de  fined  ,  i . e 


_ a. _ i. _ a. _ A.- 


(  e  [o,i] 
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It  is  further  assumed  '■hat  every  corner  point  of  3L  and  every  point  of  [■! 
at  which  there  is  a  change  in  boundary  condition  type  is  an  element 
boundary  node.  Defining 


N.(e,p)(P)  i  N.(p)U(e)(P)),  V  P  e  {P  :  P  =  P(e)(0,  £  e  [0,1] 


(8.90) 


it  follows  that 

p+1 


I  N.(e,p)(p) 
i=l  1 


=  1 


and 


N  (e,P)(p(e,p,j))  = 
i  ij 


where 

p(e  .P .  j )  r  (P^e)(c(P. j))  with  c(P.j)  =  (j-l)/p,  j  =  1.....P+1. 

The  set  of  nodes  p^e,P*j^(  j  =  l,...,p+l  consists  of  the  end  point  nodes 
of  the  element  e  and  p-1  internal  nodes. 

Let  the  set  of  all  end  point  and  internal  nodes  P  e,p’J  j  be  numbered, 
in,  say,  anticlockwise  order  around  dQ,  1  ,  . . . ,np  (see  fig.  8),  and  let  the 
coordinates  of  the  £. th  node  be  denoted  by  Thus  p^e,P’j^  =  when 


i  =  H(e  ,  j )  =  (e-1) -p  +  j . 

The  functions  and  v^ty  are  defined  in  terms  of  the  functions  N^e,P^ 


in  the  following  way: 


(u) 


Take  n^  (the  number  of  v  functions)  to  be  n-p  and  define 


1  0  • 

V"  J 


unless  P  belongs  to  an  element  having  >  as  a  node. 

(8 

when  £  is  the  jth  node  of  e  and  P  f  o . 


.91) 


Take  r>2  (the  number  of  v  functions)  to  be  n'(p+l)  and  associate 
each  v^^  function  with  an  element  and  a  node  belonging  to  the  element. 


(t)/T,.  (t) 

V  (P)  “  V(e,  j)(‘  }  ~ 


Nje,p)(P), 


unless  P  6  e, 


when  P  6  e, 


(3.9?) 


where  4'(e,j)  =  (e-l)‘(p+l)  +  j.  Note  that  the  novelty  of  defining  the 

functions  this  way  represents  an  accurate  modelling  technique  of  the 

'corner-problem'  (Fig.  9)  which  does  not  require  the  equality  of  the 

tractions  t^^  and  t^^  on  an  interelement  corner  node  (note  that 

Lachat  [44]  and  Lachat  and  Watson  [45]  assume  t.^^  =  t.^^) 

r  _  ii 

L  r  y  •  ;  ' 

It  follows  from  these  definitions  that 


(u)  . 

v  is  continuous, 

(u),D(/K  * 

Vk  (P  5  “  6kf 

r  v^u)(p)  -  l  . 

i  1 


(8.93) 


It  follows  also  that  is  zero  except  on  the  element  (or  elements)  having 

£  as  a  node,  that  v  .  x  is  zero  except  on  element  e  and  that,  when 

*  V  G  ,  J  ^ 

£'(e,j)  is  an  interelement  node,  i.e.  j  =  1 ,  or  p+1 ,  then  v  ..  is 

l  \  e  ,  j  / 

discontinuous  at  the  node  i'(e,j). 

With  these  choices  of  v^U^  and  v^ r ^  the  total  number  of  unknowns  u 

l 

and  t^,  is  n p+n(p+l).  The  set  of  points  at  which  the  integral  equation 
is  to  be  satisfied  is  taken  as  the  m  =  np  boundary  nodes  P^\...,P^m^ 

(c.f.  eq.  (8.87)).  By  imposing  the  boundary  conditions  at  these  nodes  and, 
where  appropriate,  continuity  conditions  on  t  at  interelement  nodes,  as 


many  (linear)  equations  as  there  are  unknowns  are  obtained. 


A . 1 . 3  Boundary  and  continuity  conditions 

The  interelement  nodes  are  chosen  so  that  any  point  on  the  boundary 
at  which  the  boundary  condition  type  changes  is  an  interelement  node.  At 
every  node  internal  to  an  element  it  is  required  that  the  boundary  condition 
be  satisfied.  At  every  interelement  node  the  boundary  conditions  associated 
with  each  of  the  elements  having  the  interelement  node  in  common  are  to  be 
satisfied.  Thus  there  is  associated  with  each  internal  node  two  unknowns 
(one  u^  and  one  t^,  coefficient)  and  two  equations  -  one  associated  with 
the  integral  equation  and  one  with  the  boundary  condition.  At  an  inter¬ 
element  node  there  are  associated  three  unknowns  (one  u  and  two  t  , 

Z  C 1 

coefficients)  and,  in  general,  three  equations  -  one  associated  with  the 
integral  equation  and  two  with  the  boundary  conditions. 

The  exceptional  case  arises  when  the  node  is  an  interelement  node 
between  two  nodes  which  are  both  of  Dirichlet  type,  when  the  imposition  of 
the  boundary  conditions  provides  only  one  equation.  If  the  boundary  is 
smooth  at  such  an  interelement  node  we  require,  in  addition,  that  t  be 
continuous  (and  hence  the  two  t^,  coefficients  associated  with  the  node  be 
equal)  thus  providing  one  additional  equation.  If  the  interelement  node  is 
a  corner  node,  then  from  the  directional  derivatives  of  u  along  the  two 
tangents  at  the  node  the  gradient  and  hence  normal  derivatives  of  u  can  be 
computed.  In  this  case,  then,  the  boundary  conditions  determine  the  values 
of  the  three  unknowns  associated  with  the  node  and  we  simply  exclude  from 
our  final  set  of  equations  the  integral  equation  associated  with  the  node. 

In  summary,  there  are  sets  of  equations  associated  with 

(a)  the  integral  equation, 

(b)  the  boundary  conditions,  and 

(c)  the  continuity  of  the  normal  derivative, 

and  there  are  as  many  of  these  equations  as  there  arc  unknowns. 


The  contribution  from  element  e  to  M 


k£(e ,  j)  18 

/  U(p(k),  p(e)(0)N(P)(0  — -  (O  d C 

0  J 


and,  since  v£U\p^)  =  6,  ,  the  contribution  from  element  e  to  T  ,  .. 

k  k  £  k  £  ( e ,  j ) 

for  £(e,j)  #  k,  is 

/  T(P(k),  P(e)(C))N (?)(U  £7  (£)d£, 

0  J  ai> 

where 


ds  dP 


(e) 


d£  dC 


<0 


In  both  cases  the  integrals  can  be  evaluated  accurately  using  appropriate 
gaussian  quadrature  schemes.  An  immediate  consequence  of  (8.93)  is  that 
with  ■  const,  and  t^,  =0,  u  and  £  provide  an  exact  solution  of  the 
integral  equation  (8.84).  Thus  from  (8.87)  it  follows  that 


l  T«*° 


and  thus  T^  can  be  computed  as 


Z  T 
£j*k 


k£ ' 


The  boundary  and  continuity  conditions  are  used  to  eliminate  variables 
as  the  assembly  progresses  so  that  when  assembly  is  completed,  a  system 
of  m  equations  for  the  m  unknowns  is  obtained.  When  these  have  been  solved 
the  remaining  unknowns  are  computed. 


44 


4 . 2  Treatment  of  singularities 

4.2.1  Singularities  in  a  model  problem 

In  [43]  the  following  model  problem  is  considered: 
Find  u  such  that 

V2u  -  0 


in  the  14«14  square  region  EFGH,  in  which  there  is  a  slit  (crack)  OA 
joining  the  centre  of  the  square  0  (the  crack  tip)  to  the  midside  A  of  EH 
(see  fig.  10)  with 


u  =  0  on  HA 

u  =  1000  on  AF. , 

and  3u/3n  =0  on  the  remaining  boundary. 

By  symmetry  u-500  is  antisymmetric  about  the  line  AB  and  in  the  rectangular 
region  EFBOA  u  is  the  solution  of  the  'model  problem' 


(Fig.  10  here) 


V2u  =  0 
u  =  1000 
u  =  500 
3u/3n  =  0 


in  EFBOA 
on  AE 
on  BO 

on  OA,  EF  and  FB 


Using  the  standard  separation  of  variables  technique,  the  boundary 
conditions  on  OA,  the  boundedness  and  anti-symmetry  of  u  -  500,  one 
obtains  (see  e.g.  [46]) 


u(P) 


u(0) 


00 


+  ?; 
j=o 


ft  .  f  •  ( r  ,0 ) 
J  J 


in  the  neighbourhood  of  0,  and,  hence,  for p on  OB, 


00 


l  a  . 

j“0  J 


1  ?fi 

r  36  (r,0)’ 


( R . 04  ) 


( R .05) 


where 


f.(r,6)  i  r(2j  +  D/2  s in{  (--^)0) 

and  (r,G)  are  the  polar  coordinates  of  P. 

In  problems  of  thin  kind  the  coefficient ,  ,  oi  the  leading  term, 

is  of  special  interest,  e.g.  in  the  crack  problem  the  corresponding 
coefficient  is  proportional  to  the  stress  intensity  factor  (SIF).  Note 
that,  for  the  model  problem  we  have 


lim  r  "1[u(r,n)  -  u(0,0)], 
r-0 


lim  2r  —  (r,0). 

„  9n 
r-*0 


(«.%) 


(8.97) 


The  boundary  element  method  as  described  in  section  (4.1)  is  now 
applied  to  this  model  problem.  Let  the  boundary  of  the  region,  EFBOA, 
be  subdivided  into  n  elements  and  let  the  elements  be  numbered 
anticlockwise  from  0  as  shown  in  Fig.  11a. 

.  .  -i/2 

Since  on  OB,  in  the  neighbourhood  of  0,  t*  behaves  like  r  , 

while  7  is  bounded,  we  cannot  expect  to  obtain  an  accurate  representation 

1/2 

of  t*  near  0.  On  OA,  m  the  neighbourhood  of  0,  u*  behaves  like  r 

(1  2) 

while  on  element  1,  if  the  internal  node  P  ’  is  chosen  to  be  mid-way 
between  the  element  boundary  nodes,  and  p  is  taken  to  be  2,  u  will  be  a 
quadratic  function  of  r  and  hence  cannot  accurately  represent  u* . 

Thus  we  would  not  expect  the  method  of  section  4.1  to  work  well  for 
this  problem  and  in  particular  we  would  expect  poor  results  for  the 
computed  values  of  a^.  In  this  section  we  describe  various  ways  of 
overcoming  these  difficulties  and  in  section  4.3  we  attempt  to  compare 
their  efficacy.  Some  of  the  techniques  can  be  extended  in  a  natural  wav 


to  deal  with  singularities  in  more  general  problems  for  which  the  form  of 
the  singularity  is  known.  While  techniques  of  the  kind  described  here  have 


4', 


been  used  in  standard  finite  element  method  calculations  of  stress 
intensity  factors  (S1F)  for  crack  problems,  the  use  of  some  of  then 
in  the  BIE  context  appears  to  be  new.  In  the  standard  FE  method 
calculations  of  SIF's  basically  three  types  of  methods  have  been  used: 

(1)  Those  methods  which  use  special  functions  having  the  correct 
singular  behaviour  in  the  neighbourhood  of  the  singularity  (see  sections 
4.2.2,  4.2.3).  For  the  BIE  such  special  functions  can  be  introduced 


by  modifying  the  definitions  of  v^  and  v^  of  section  4.1  to  include 


(u) 


functions  having  the  correct  singular  behaviour.  For  v  ,  an  elegant 


alternative  is  to  choose  element  internal  nodes,  P^e’^,  so  that  v^U^ 


1/2 


has  the  r  behaviour  of  u. 

(2)  Those  methods  which  transform  the  problem  into  one  in  which 
singularities  are  not  present.  Two  such  methods  are: 

(a)  subtraction  of  singularity  techniques  (see  section  4.2.4),  and 

(b)  conformal  mapping  techniques.  This  technique  has  been  used 
successfully  for  problems  in  two  dimensions  (see  e.g.  [47]).  However, 
since  it  cannot  be  extended  to  deal  with  problems  in  three  dimensions 
it  is  excluded  from  our  discussion. 

(3)  Those  methods  which  do  not  use  special  functions  but  employ  fine 
meshes  near  to  the  singularity. 

From  a  computed  solution  obtained  by  whatever  method,  one  can  compute 


approximations  to  a q  by  a  variety  of  schemes  (see  section  4.3).  If  the 


solution  is  calculated  by  a  given  method,  for  a  succession  of  grid 
refinements,  then  'extrapol ation  to  zero  grid  size'  of  the  values  of  a 


0 


computed  by  one  such  scheme,  may  yield  an  excellent  approximation  to 


o’ 


see  e.g.  Morley  [48]  where  extrapolation  leads  to  a  considerably  improved 
approximation . 

Extrapolation  techniques  will  give  good  results  when  (i)  the  error, 


e(h),  in  the  computed  value  of  behaves  asymptotical ly  like  some  power 


of  h,  where  h  is  a  measure  of  fineness  of  the  grid,  and  (ii)  the  computations 


47 


are  carried  out  for  values  of  li  for  which  e(h)  is  well  approximated  by 
its  asymptotic  form.  When  condition  (i)  holds  one  would  expect  that  grids 
for  which  e(h)  has  this  property  will  be  relatively  much  coarser  for  those 
methods  for  solving  boundary  value  problems  in  which  special  provision  is 
made  for  treating  singularities,  than  in  those  methods  in  which  no  such 
provision  is  made. 


4.2.2  Treatment  of  singularities  in  u*  and  t*  using  special  functions 

-1/2 

(A ) .  A  treatment  of  the  r _ singularity  in  t*  using  special  functions 

-1/2 

Suppose  that  we  wish  to  represent  the  r  behaviour  over  an  interval 

OB'  where  B'  lies  on  the  interval  OB.  The  element  subdivision  is  chosen 

so  that  B'  is  an  interelement  node  and,  for  each  element  belonging  to  OB', 

(e  2) 

the  internal  element  nodes  P  ’  are  taken  to  be  mid-way  between  the 

(el)  (e  3) 

element  boundary  nodes  P  ’  and  P  ’  so  that 

P^C^(f)  =  p^e’^  +  £(p^e’^  -  p^e’^) 

(e) 

Let  the  length  of  the  element  e  be  £  and  let  the  radial  coordinates  of 

n(e,l)  T^e),.  ,  ,(e)  .  (e),.,.  .  .  , 

B  ,  P  and  P  (C)  be  r  ,,  d  and  r  (£),  respectively  (see 

D 

fig.  11(b)).  Thus  for  C  6  [0,1] 


r(e,<5)  -  l(e)(»!e)  -  {). 


where 


A(e)  =  d(eVe) 


We  define  v^i(e  j)  as  in  section  4.1,2  except  that,  for  e  €  OB',  we 
replace  in  the  definition  (8.89)  by  where 


N(p.’e)(C)  =  (r./r(e)U))1/2N(p)U). 

IB  1 


I 
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(t) 


Thus,  this  definition  of  v„,,  ..  differs  from  that  given  in  section  4.1.2 

l  ( e  ,  j  ) 


only  when  e  6  OB'.  Note  that  the  i nterelemen t  continuity  between 
sxngul ar/continua tion/ordinary  boundary  elements  is  preserved.  (Fig.  11(c)) 
It  follows  at  once  from  this  definition  that ,  for  every  set  of 


coefficients  t.,,  .»  for  which  the  continuity  conditions  of  section  4.1.3 

*■  (e,  j) 


”1/2 

are  satisfied,  on  OB',  t(P)  has  the  form  P^(r)r  ,  where  (r,0)  are  the 


polar  coordinates  of  the  point  P,  P^(r)  is  continuous  and  is  a  pth  degree 


polynomial  on  each  element  e;  conversely,  to  every  such  function  P^(r) 


there  corresponds  a  set  of  coefficients  t .  ,  .  .  .  such  that  t  has  the  form 

t'(e.j) 


Pj(r)r  Thus,  using  these  singular  v^^  functions  we  can  accurately 


represent  the  sum  of  the  first  (p+l)  terms  in  the  expansion  (8.93)  for 


* 

t  on  OB' . 


/  1/2  * 

(B)  A  treatment  of  the  r _ singularity  in  u  using  special  functions 


(1  2) 

Let  the  internal  node  P  ’  of  element  1  be  mid-way  between  the 


element  boundary  nodes  and  P^’^  and  let  the  functions 


be  defined  by 


S(f<0  -  u(?,/{)1/y?>«:>. 


i  =  ; , . . . ,p+i 


N(^(0  =  1  -  ?l  N*?*(0 

i  =  2  1 


We  define  as  in  section  V.hAbut  with  N^?^(0  replaced  by  N^?^(£) 


(u) 


when  e  *  1 .  Note  that  in  this  definition  the  v  retain  the  properties 


(8.93).  It  follows  from  this  definition  that,  on  element  1,  u  has  the  form 


1/2, 


u(P)  =  u(0)  +  r  Q^(r) 


where  r  is  the  radial  coordinate  of  the  point  P  and  Qj(r)  is  a  (p-l)th 


(u) 


degree  polynomial  in  r.  Thus,  using  the  functions  v  ,  defined  in  this 


way,  u  can  accurately  represent  the  sum  of  the  first  p  terms  in  the 


_ _ _ _ «-  •  a-*  -  I 


3|( 

expansion  (8.94)  for  u  ,  on  the  element  1.  Unfortunately,  there  does  not 

seem  to  be  any  simple  way  of  introducing  special  functions,  N,  for 

1/2 

elements  2,3 .  to  reproduce  the  r  behaviour  over  those  elements 

while  retaining  the  property  (8.93). 


1/2 

4.2.3  Treatment  of  the  r _ singularity  in  u  using  the  element 

parametrization 

1/2 

An  alternative  way  of  treating  the  r  singularity  in  u  which 
does  not  involve  the  introduction  of  special  v^U^  functions  was  suggested 
by  a  device  originally  proposed  independently  by  Barsoum  [49]  and  Henshell 
and  Shaw  [50]  and  extended  by  Lynn  and  Ingraffea  [51]  in  the  context  of 
the  finite  element  method.  J  j  I  ?  ,  ’ 

Consider  a  typical  element  e  of  length  2h  on  OA  as  shown  in  fig.  12. 

(e  2)  (el)  (e)  (el) 

Let  the  distance  of  P  ’  from  P  *  be  p  'h  and  the  distance  of  P  ’ 

(e) 

from  the  point  0  be  q  h.  Let  r  be  the  polar  coordinate  of  the  point 

(e)  (e  2) 

P  (O  and  let  the  position  of  the  internal  node  P  ’  be  chosen  so  that 


(e)  _  1  n  (e)  /(e).  (e)  ,,, 

p  =  r  (1-q  +  (q  +  2))  , 


(8.98) 


then 


£  =  Q1  (q^G^){l  -  /r/(q^h)  )  , 


where 


QL(q)  =  l/{  1  -  </(q  +  2)/q) 


^  .  .  1  /2 

and  hence  IT  will  be  a  pth  degree  polynomial  in  r  on  the  element  e. 

Thus  if  the  position  of  the  internal  nodes  of  the  elements  1,2,..., 

belonging  to  OA  are  chosen  according  to  (8.98),  then  "u  can  accurately 

represent  any  continuous  function  which  is  piecewise  a  pth  degree'  polynomial 

in  r  on  those  elements.  In  particular,  taking  p  =  1,  u  can  accurately 


represent  any  function  of  the  form 


5o 


u(o) 


Thus,  using  this  scheme,  with  p  =  1 ,  we  can  reproduce  the  first  term 
in  the  expansion  (8.94)  for  u  near  0. 

Unfortunately,  to  reproduce  the  first  q  terms  in  the  expansion,  we 
must  take  p  =  2q  -  1 . 

In  the  boundary  element  method  described  in  section  4.1  the  boundary 

of  the  region  is  approximated  by  a  curve  which  is  defined  parametrically 

in  terms  of  quadratic  Lagrangian  interpolation  polynomials.  An  obvious 

generalization  of  this  specification  of  the  boundary  is  obtained  by  taking 

m  >  1  internal  nodes  per  element  and  using  (m+l)th  order  Lagrangian 

polynomials  in  the  parametric  definition  of  the  boundary.  Pu  ,  Hussain  and 

Lorensen  [52]  have  extended  the  work  of  [49]  and  [50]  to  the  case  m  =  2. 

They  give  formulae  for  the  position  of  the  two  internal  nodes  on  element  1 

1/2 

such  that  £  is  linear  in  r  on  the  element.  Bernal  and  Xanthis  [53]  have 

derived  formulae  for  the  position  of  the  two  internal  nodes  on  any  element 

1/2 

belonging  to  OA,  such  that  £  is  linear  in  r  on  all  elements  e  6  OA . 

For  the  case  m  =  2,  with  this  choice  of  internal  nodes,  v^U\  as  defined 

1/2 

by  (8.91),  has  the  desired  r  behaviour.  However,  results  are  presented 

-1/2 

for  the  case  r  =2,  p  =  2  only.  In  summary,  we  can  introduce  the  r 

behaviour  of  t*  into  t  by  using  the  modified  Lagrangian  interpolation 

functions,  N^!e\  and  the  r^^  behaviour  of  u*  into  u  by  using  the  modified 

'(p) 

Lagrangian  functions,  N  5.  ,  or,  for  the  case  p  =  2 ,  by  choosing  the  internal 
(e  2  ) 

element  nodes,  P  ’  ,  (for  e  -  1,2,...,  and  e  6  OA)  according  to  (8.98). 

An  alternative  way  of  treating  the  singularity  is  to  transform  the  problem, 
using  a  subtraction  of  singularity  technique,  so  that  the  solution  oi  the 
transformed  problem,  to  be  solved  by  the  BIE  method,  does  not  possess 
singularities. 


4.2.4  Subtraction  of  smculari  ties 


In  the  neighbourhood  of  0 


u(r,0)  =  u(0)  +  Z  a.f.(r,9), 

j=0  J  J 

where 

f.(r,0)  =  r(2J+1)/2  s in{y( 2 j+I )9[ 

Thus  we  might  reasonably  expect  that,  near  0,  u  could  be  approximated 
closely  by  retaining  only  a  few  (say,  q+1)  terms  of  the  sum,  i.e.  that 


u=u(o)+  I  a.f.  ,  nearO. 
j=0  J  J 


Let  us  define 


i/q)  =  u  -  Z  a.f.  , 

j-0  J  J 


(8.99 


then  , 

(i)  if  were  known  we  could  compute  u  at  once  from  (8.99), 

(ii)  since  the  f^  are  harmonic  functions,  is  the  solution  of  the 

boundary  value  problem: 

find  such  that  =0  in  EFBOA , 


1000  -  Z  a.f. 
j=0  J  J 


500  -  Z  a.f. 

j-0  J  J 


on  AE  , 


on  BO , 


q  9  f 
-  T.  a.  — 


on  OA ,  FF  and  FB, 


(iii)  whereas  u  has  d i scon t inuous  derivatives  at  0, 


1  q; 


ha  s 


continuous  derivatives  of  order  q. 

q  )  Vii 

Any  numerical  scheme  for  computing  an  approximation  ijj  M  f 

•  „(q) 

provides  an  approximation  u  to  u: 


a  .  f  • 

J  J 


and  by  virtue  of  (iii)  we  might  expect  that  an  algorithm  (e.g.  finite 
difference,  finite  element,  BIE  method)  in  which  no  special  steps  are 
taken  to  treat  singularities,  would  consequently  give  a  reasonable 
approximation  0(q)  to  the  solution  t  ^  ^  of  and  hence  that 

de f ined  by  (8.99),  would  give  a  reasonable  approximation  to  u. 

Of  course,  to  use  such  an  approach  we  would  have  to  know  a  ,  ...,a  , 

0  q 

or  at  least  be  able  to  determine  sufficiently  good  approximations  to  these 
constants.  Thus  the  equations  derived  from  whatever  scheme  we  use 
(FD,  FE ,  etc.)  must  be  supplemented  by  additional  equations  such  that  the 
total  system  of  equations  determines  simultaneously  an  approximation  0^^ 

to  and  approximations  to  a.  (j  =  0 . q).  Since,  as  we  have 

al ready  noted , 

q 

u  -  u(0)  +  E  a.f.,  near  0, 

j=0  J  3 

we  might  take,  in  the  case  of  finite  difference  or  finite  element  schemes, 
the  additional  equations  to  be,  for  example,  those  obtained  from  the 
requi remont 

~<q)  =  u(0)  , 


at  (q  +  1)  nodes  near  0,  not  belonging  to  that  part  (  rif;  )  of  the  b^undirv 
On)  on  which  u  is  prescribed. 
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In  the  case  of  the  BIE  method  described  in  section  4.1  with  the 
approximations , 


T(q)‘  = 


y  .  (U) 

zv*  > 


-(q) *  _  y  ( t ) 
x  -  ZxrV 


(q)  *  (q)*_  (q) 

to  if  and  X  =  3C'  M  /3n,  respectively,  we  might  take  as  additional 

equations  those  obtained  from  the  requirement  that  be  continuous 

at  0,  =  u(0)  at  Tj  points  of  3fi  near  0,  not  on  3fl  ,  and  =  0 

at  r2  points  of  3Q  near  0  and  not  on  that  part,  3Q2,  of  the  boundary  on 
which  3u/3n  is  prescribed,  where  -  q. 

Thus,  with  elements  numbered  as  in  fig.  11,  and  with  p  =  2 ,  we 
might,  for  example,  take 


X(n,3)  =  °  if  q  =  °* 


and  additionally 


x(n  2)  =  0  ’  ^2  =  u('0^  if  q  =  2  > 


and  additionally 


(n ,  1 ) 


=  0,  i^3  =  u(0)  if  q  = 


The  subtraction  method  outlined  above  can  easily  be  extended  to  deal 
with  any  problem  for  which  the  form  of  the  singularity  is  known  (e.g. 
any  reentrant  comer  problem).  Symm  [54]  was,  perhaps,  the  first  to  use 
this  technique  in  the  context  of  integral  equations  (his  approximating 
functions  v^U^  and  v^^  were  step  functions).  See  also  Papamichael  and  Symm 
[55]  where  the  work  of  [54]  is  amplified.  Xanthis  et  al  [43]  hnv  ii'o1  th 
subtraction  of  singularity  technique  with  both  the  finite  difference  r>  th.-l 
and  the  BIE  method  described  in  section  4.1.  In  both  cases  the  ri-tk-H 


worked  well.  Their  BIE  calculations  arc  reported  in  the  following  =■'<-? 


5'. 


The  program  for  computing  the  coefficients  if’  and  >^ ,  and  the  coefficients 
a ^ ,  j  =  0,...,q,  by  the  BIE  method,  was  obtained  by  a  relatively  trivial 
modification  of  the  program  described  briefly  at  the  end  of  section  4.. 4. 


A . 3  Numerical  results  for  the  model  problem 

In  [A3],  the  basic  BIE  method  described  in  section  A.l  and  modifications 
of  this  method,  using  the  techniques  described  in  section  A. 2  have  been 
used  to  compute  solutions  of  the  model  problem  given  in  Fig.  10,  of 
section  A . 2  . 1 . 

Each  of  the  five  boundary  segments  OA,  AE ,  EF ,  FB  and  BC  were 
subdivided  into  N  elements  of  equal  length,  and  for  N  =  2,  3,  A  and  5  the 
solution  of  the  model  problem  was  computed  by  each  of  the  following  modifi¬ 
cations  of  the  basic  method  described  in  section  A.l: 

Method  1.  No  modification. 


Method  2.  v^U^  functions  with  r^“  behaviour  on  element  1  (section  A. 2) 


(t)  .  .  -1/2 

v  functions  with  r  behaviour  on  OB  (sections  A. 2. 2). 


Method  3,  nodes  on  OA  chosen  so  that  v^U^  has  r^2  behaviour  on  OA 


(t) 


(section  A. 3);  v  as  in  Method  2. 

Method  A ,  subtraction  of  singularity  technique  (section  A. 2. A). 
For  this  method  the  computations  were  carried  out  for  q  =  0,  2  and  A. 

For  each  computed  solution,  the  values  of  u  at  the  points  (r,-n), 
r  =  1/A  and  r  =  3/A,  were  calculated. 


A. 3.1  Computed  u  values 

These  are  given  in  tables  la  and  lb.  From  these  tables  it  is  evident 
that,  for  a  given  value  of  N,  method  A  with  q  =  A  gives  the  most  accurate 
results;  method  2  and  method  A  with  q  =  2  give  results  which  are  nearly  as 
accurate  as  those  obtained  from  method  A  with  q  =  A.  However,  this 
conclusion  does  not  imply  thatw^,  the  parameter  which  is  usually  of  most 
importance  in  the  applications,  would  also  be  given  most  accurately  by 


me  thod  A . 


Table  la 

Values  of  u  at  r  =  1/4;  'exact*  value  =  576.41  (541 


ri 

Me  thod 

1 

2 

3 

4 

N 

q  =  0  q  =  2 

534.17 

576.22 

574.79 

576.51 

376.48 

576.44 

540.69 

576.31 

575.69 

576.48 

576.42 

576.41 

545.96 

576.36 

576.06 

576.47 

576.41 

576.41 

550.38 

576.38 

576.25 

576.42 

576.41 

576.41 

Table  lb 

Values  of  u  at  r  -  3/4:  'exact*  value  =  634.45  [54] 


Method 

1 

2 

3 

4 

N 

q  =  0 

596.96 

634.29 

634.15 

634.68 

634.56 

634.50 

611.49 

634.40 

634.77 

634.50 

634.46 

634.45 

621.29 

634.44 

634.88 

634.43 

634.45 

634.45 

627.70 

634.45 

634.85 

634.45 

634 .45 

634.45 

A . 3 . 2  Calculation  of  an  approximation  t_o 

For  Method  A,  above,  we  can  take  •  By  whatever  method  u 

and  t  are  computed,  one  can  use  the  formulae  obtained  by  substituting 
u  and  If  into  the  r.h.s.  of  (8.96)  and  (8.97),  respectively,  to  obtain 
approximations  to  a^: 


=  a^^(u)  -  lim  r  ^2[u(r,ir)  -  u(0)] 
r-*0 


(8.96) 


and 


aQ  =  c/^(t)  =  lim  2r^“t(r,0). 
r-HD 


(8.97) 


Note  that  if  u  and  t  are  obtained  from 
( i)  Method  A ,  then 


(u),-.  (t),~. 

a  0  (u)  -  a  0  (t)  -  V 


(ii)  Method  1,  then 


/  *  ( U  )  (t),r\  o 

a  Q  (u)  =  a  Q  (t)  =  0. 


Nevertheless,  in  case  (ii),  one  can  obtain  useful  approximations  to 
from  u  and  t  using  the  approximations  to  (8.96)  and  (8.97): 


a^U^(u,r  )  =  r  ^2[u(r  ,tt)  -  u(0)] 
0  u  u  u 


(8.96) 


and 


n^(t\rt)  =  2r^2f(rf  ,0) 


(8.97) 


with  ry,  rt  sufficiently  small.  Formulae  (8.96)^  and  (8.97)^  can,  of 
course,  be  used  to  compute  approximations  to  however  u  and  t  are 
calculated.  For  Method  1,  experiments  showed  that  there  exist  constants 


a,  v,  a*,  v*  such  that,  with 


57 


r  =  a£  , 
u 


=  a*?V 


(where  i  =  length  of  element  1  =  length  of  element  n  =  7/N), 

,r^)  and  a^t^(T,rt)  approximate  closely  aQ.  The  attempt  to 

determine  such  constants  was  prompted  by  the  observation  of  Schatz  and 

Wahlbin  [56]  that,  for  the  standard  finite  element  method,  a  good  choice 
/  \  /  *\ 

of  r  for  a  q  is  n  where  h  is  the  maximum  diameter  of  all  elements. 

The  fact  that  we  were  able  to  find  constants  a,  v,  a*,  v*  so  that  a0  is 
closely  approximated  by  o/^Ci^r^)  and  a  ^(t.r^)  i-s  obviously  significant 

and  itcurrently^/investigatoJ.  . _ __ _ —  Additional 

approximations  to  aQ  can  be  obtained  from  formulae  deduced  from  certain 
path-independent  (invariant)  integrals,  the  F-  and  M-integrals,  evaluated 
along  a  circuit  surrounding  the  point  where  the  singularity  occurs  (see 


section  2.2).  These  integrals  are  defined  by 


Ft  ‘  I  Vi  “S' 


M  =  f  x  P  .  n .  dS 
^  t  jt  j 


(8.100) 


(8.101) 


iere  S  is  a  plane  curve  with  normal  n  =  (n^n^),  the  energy  momenti 


tensor 


V  ‘  U£j  '  ~ .  ",l 

>  3 


12  2 

and  L,  the  Lagrangian  function,  L  =  y(u  ^  +  u  ■ 

Taking  the  x2~axis  to  be  along  AB  and  S  to  be  a  closed  path  containing 
the  point  0  -  the  crack  tip  (see  fig.  10)  -  one  obtains 


V2  ’ 


aTidg/2  . 


(P.102) 


M 


(8.103) 


In  particular,  taking  S  to  be  the  square  EFGH  and  using  the  symmetry 
properties  of  the  integrals  and  the  computed  u  and  t* values,  one  can 
find  approximations  and  M  to  F2  and  M,  and,  hence,  approximations  to 
cIq  from  (8.102)  and  (8.103),  which  are  denoted  by  a^^(u,T)  and  a^(~,t') 
respectively. 

From  the  u  and  t  values,  computed  by  the  four  BIE  variants,  Methods  1, 
2,  3  and  4,  approximations  to  are  calculated  and  the  results  given  in 
the  following  tables. 

Inspection  of  table  2b,  c,  d  and  e,  shows  that  for  the  Methods  2,  3 
and  4  one  can  obtain  very  accurate  values  for  ,  with  a  coarse  grid 
(N  =  2,3,4),  using  the  approximations  c/^(*u,t),  a^^(t)  and,  for 
Method  4,  a^.  For  Method  1  to  give  an  approximation  to  using  a  Q  (u,t) 

-(M)  s 

or  a  (u,t),  to  an  accuracy  comparable  to  that  obtained  from  Methods  2, 

3  and  4  with  a  coarse  grid  (N  :  5),  a  much  finer  grid  must  be  used.  For 
example,  N  =  10  gives  a^\u,t)  =  151.41  and  o/^fu.t)  =  151.76,  values 
which  are  not  as  accurate  as  those  obtained  by  the  other  methods  with 
N  =  5.  Nevertheless,  the  values  a  ^  (u,t)  and  a  ^  (u,t)  obtained  using 
Method  1  with  N  =  5  are  accurate  to  approximately  one  per  cent.  Perhaps 
it  is  worth  noting  that  (i)  (t)  is  obtained  as  a  simple  multiple  of 

one  of  the  unknowns,  whereas  a^\u,t)  and  c/^(u,t)  involve  the  computation 
of  line  integrals,  and  (ii)  Method  4  gives  without  any  further  compu¬ 
tations,  although  it  must  be  remembered  that  the  system  of  equations  from 
which  the  solution  is  obtained  involves  q+1  more  equations  than  the  other 
methods.  Extrapolation  to  zero  grid  size,  for  a^^(u),  was  carried  out  by 
assuming  that 

+  A(l/N)P  , 


the  constants  a^,  A  and  8  in  this  formula  being  determined  from  the  values 
of  an(u)  for  three  successive  values  of  N.  For  method  2,  using  N  =  2,3.4 
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and  N  =  3,4,5,  the  computed  values  of  were  151.78  and  151.64 
respectively.  For  Method  3  the  corresponding  results  were  151.58  and 
151.64,  respectively.  These  results  demonstrate  clearly  the  benefits  of 
extrapolation. 

4.4  Concluding  remarks 

An  improved  implementation  of  the  BIE  method  has  been  discussed, 
which  provides  a  special  treatment  of  discontinuities  associated  with 
corners  and  points  at  which  there  is  a  discontinuity  in  the  boundary 
condition.  Three  variants  of  this  basic  algorithm  have  been  introduced 
with  modifications  capable  of  accommodating  boundary  singularities. 

It  must  be  noted  that,  although  one  boundary  singularity  has  been 
considered,  any  number  of  such  singularities  may  be  treated  simultaneously. 
The  same  treatment  of  the  square-root  singularity  that  has  been  described 
in  the  context  of  harmonic  problems  is  directly  applicable  (component 
wise)  to  similar  situations  involving  Navier's  equation  of  elasticity. 

An  attractive  feature  in  the  three  methods  described  in  section  4  is  that 
the  stress  intensity  factor  appears  as  one  of  the  unknowns  in  the 
computation  or  a  simple  linear  combination  of  these  unknowns,  thus  avoiding 
the  uncertainties  associated  with  schemes  which  involve  plotting  and 
extrapolation  techniq  es. 

In  addition  to  the  improved  implementation  the  results  of  [43]  show 
how  the  invariant  integrals  of  section  2.2  can  be  used  both  directly  and 
as  a  check  on  other  methods  of  stress  intensity  factor  determination. 

These  methods  have  also  been  applied  to  certain  crack  problems  in  elastic 
media  with  spatially  varying  elastic  properties  [57]. 


Method  4  Method  3  Method  2  Method 


Table  2 


Computed  values  by  three  BIF  variants.  For  Method  4  higher  order 

coefficients  for  a.  are  also  presented.  'Exact'  value  of  =  131.63 
J  0 


N 

2 

3 

4 

5 

(F),  ^ 

~  a  0 

155.52 

154 .22 

153.57 

153.18 

Method 

Q 

O  3! 

v-/ 

c 

r* 

152.29 

152 .07 

151.96 

151.89 

(a) 

N 

a  0  (^) 

wsm 

mm 

(M)r  77 

a  0  (u,t) 

2 

151.13 

151.74 

151.66 

151.57 

3 

151.35 

151.65 

151.63 

151.62 

(b) 

4 

151.46 

151.63 

151.63 

151.62 

5 

151.51 

151.63 

151 .63 

151.62 

N 

(u) 

a  Q  (u) 

(t),-. 
a  0  (t) 

( F) . - 
a  0  (u,t) 

(  M )  /  -  -  -. 

a  0  (u  ,t ) 

2 

142.32 

151.72 

151.66 

151.56 

3 

145.58 

151.64 

151.63 

151.61 

(c) 

4 

147.17 

151.63 

151.63 

151.62 

!  5 

148.11 

151.63 

151.63 

151.62 

N 

q 

o  q 

ao  ai 

°2 

2 

150 

.56 

151.71  4.68 

0.13931 

,  3 

0  150 

.86 

151.63  4.71 

0.13838 

(d) 

:  4 

0  151 

.05  2 

151.63  4.72 

0.13690 

!  3 

151 

.17 

151.63  4.73 

0.13611 

N 

q  a 

o  ai 

a  „ 

2  3 

a4 

-3* 

2 

151.70 

4.73 

0.13069 

0.00866 

0.00027 

•o 

3 

151.63 

4.73 

0.13296 

0.00883 

0.00023 

o 

X 

4 

4  151.63 

4.73 

0.13291 

0.00886 

0.00023 

V 

V 

5 

151.63 

4.73 

0.13295 

0.00887 

0.00022 
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1:  Definition  of  crack  modelling  terms  used 

in  deriving  equation  (  .16)  (After  Cruse  [3,4,5]) 

2:  Typical  crack  tip  co-ordinate  system 

3:  Fracture  modes.  3(a)  opening  mode;  3(b)  sliding 

mode;  3(c)  anti-plane  shearing  mode 

4:  Contours  for  integrals  around  crack  tip 

5:  A  crack  positioned  symmetrically  in  a  body 

6:  A  crack  positioned  unsvmmetrical ly  in  a  body 

7;  A  Somigliana  dislocation 

8:  Numbering  of  elements  and  nodes  :  case  p  =  4 

9:  An  interelement  comer  node 

10:  Square  region  with  a  slit 

11:  (a)  Element  numbering  for  the  model  problem 

(b)  Notation  used  in  modelling  t  ahead  of  the  crack 

(c)  Illustration  of  singular  and  first  (compatible) 
continuation  Lagrangian  'square-root'  elements 

12:  Notation  used  in  modelling  u  on  the  crack 
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